Geochemical modeling of a sustained shallow aquifer CO2 leakage field study and implications for leakage and site monitoring

International Journal of Greenhouse Gas Control 37 (2015) 127–141

Contents lists available at ScienceDirect

International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc

Geochemical modeling of a sustained shallow aquifer CO2 leakage field study and implications for leakage and site monitoring Aaron G Cahill a,∗ , Rasmus Jakobsen b a b

G360 Centre for Applied Groundwater Research, School of Engineering, University of Guelph, Ontario, Canada Geological Survey of Denmark and Greenland, Øster Voldgade 10, 1350 Copenhagen K, Denmark

a r t i c l e

i n f o

Article history: Received 2 April 2014 Received in revised form 9 March 2015 Accepted 13 March 2015 Keywords: Carbon capture and storage Carbon dioxide Leakage Groundwater chemistry Trace metals Reactive transport modeling

a b s t r a c t A geochemical numerical modeling study was conducted to constrain processes occurring in field and laboratory experiments, simulating CO2 leakage from geological storage on shallow potable aquifers. A leak was previously physically simulated in a shallow potable aquifer at Vrøgum plantation, Western Denmark by injection of 1600 kg of gas phase CO2 over 72 days. Here, a 1-dimensional reactive transport model was constructed based on field and laboratory results and subsequently used to explore the contributions of various geochemical processes to explain observed results from the carbonate free system. Finite gibbsite derived Al3+ driven cation exchange is able to explain the majority of water chemistry change observed at Vrøgum including: a pulse like effect showing a fast peak and return toward background levels for alkalinity and dissolved ion concentrations; and increasing and persistent acidification via buffering exhaustion. Model processes were supported further by simulation of a batch experiment conducted on the Vrøgum glacial sand, employing the same processes and sediment parameters. The fitted reactive transport model was subsequently used to extend predictions and explore various scenarios. Extended predictions suggest the pulse of elevated ions travels with advective flow succeeded by a zone of increasing acidification. Model runs at higher PCO2 (implying greater depths) suggest amplification of effects, i.e., greater peaks and more rapid and severe acidification. Calcite limits acidification, however, induces additional Ca driven ion exchange giving rise to more significant chemistry change. Although a site specific model, results have significant implications for risks posed to water resources from CCS leakage and implementation of MMV programs. © 2015 Elsevier Ltd. All rights reserved.

Introduction Assessment of safety and feasibility of carbon capture and storage (CCS) in deep saline aquifers in order to mitigate climate change is ongoing and the focus of many research teams throughout the world. For public acceptance and global implementation, comprehensive understanding of leakage from CCS is required (Boyd et al., 2013). Assessment of leakage requires consideration of all physical and chemical processes occurring in the reservoir, cap rock, overlying fresh water formations and ultimately atmospheric discharge. Over the past 10 years all aspects of leakage have been investigated and one key focus has been on the effects elicited by CO2 on water quality in shallow potable aquifers overlying storage reservoirs. The shallow aquifer focus can be broadly described as having two main aims; (1) to determine how deleterious a

∗ Corresponding author. Tel.: +1 519 8244120x56269. E-mail address: [email protected] (A.G. Cahill). http://dx.doi.org/10.1016/j.ijggc.2015.03.011 1750-5836/© 2015 Elsevier Ltd. All rights reserved.

leak may be to groundwater resources and human health (Siirila et al., 2012) and (2) how best leakage can be detected geochemically (for monitoring purposes) (Klusman, 2011). In decreasing number; numerical modeling, field studies and laboratory studies have all recently been used to characterize the likely effects of CO2 leakage on shallow aquifers (Lemieux 2011; Lions et al., 2014). Consequently the effects of CO2 on water chemistry and sediment are becoming increasingly understood, i.e., reduction in pH, increases in major and minor ions, dissolution/precipitation of reactive minerals and sorption/desorption processes. Modeling studies for the most part have taken a generic approach (Jaffe and Wang, 2003; Carroll et al., 2009; Zheng et al., 2009, 2013; Wilkin and Digiulio, 2010; Humez et al., 2011; Navarre-Sitchler et al., 2013) while field and laboratory studies have explored site specific effects on aquifers of varying mineralogical and hydrogeological composition (Assayag et al., 2009; Kharaka et al., 2010; Lu et al., 2010; Cahill et al., 2013; Varadharajan et al., 2013; Cahill et al., 2014; Humez et al., 2014). Historically the number of published numerical modeling studies has been high (greater than 15) in relation

128

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

to laboratory and field studies (less than 10). Recently this imbalance appears to be equalizing as several new field studies have been reported exploring the effects of CO2 on a variety of geologic settings and temporal and spatial scales (Peter et al., 2012; Cahill and Jakobsen, 2013; Trautz et al., 2013; Cahill et al., 2014). The main challenge associated with this work is determining applicability to a real leakage scenario. Generic modeling studies lack site specific data with which to validate their predictions, laboratory experiments are not reflective of in situ conditions (i.e., flow and/or sediment:water ratios are not considered) and field studies are site specific and tend to be relatively limited in scale (spatially and temporally). For this reason an integrated approach is needed whereby field and laboratory experiments must be conducted then simulated numerically to test understanding, assess model capabilities and quantify processes before long term predictions can be made. There are currently several studies which take this approach. For example (Zheng et al., 2012) used the TOUGHREACT V2 (Xu et al., 2011) code to model results obtained from the first CO2 contamination field experiment conducted (Kharaka et al., 2010). Batch style geochemical simulations were employed disregarding flow due to insufficient data and uncertainties in flow regimes. Three models were utilized implementing combinations of key geochemical processes (dissolution/precipitation, ion exchange and surface complexation). The results showed good approximation to field data for most parameters with calcite derived Ca2+ driven cation exchange as the likely cause of changes in water chemistry with surface complexation processes being negligible. Another study (Viswanathan et al., 2012) used a combination of mineral characterization techniques, laboratory batch experiments and TOUGHREACT (both batch and reactive transport set up) to assess exchange, sorption and dissolution/precipitation processes controlling water chemistry at a natural analog for CO2 leakage in Chimayo, New Mexico (Keating et al., 2010). Initial increases in arsenic were observed followed by steady consistent decreases which were successfully simulated by including precipitation of clay minerals (illite and kaolinite) and sorption. Another modeling study also used this natural analog site to inform model development, further exploring processes taking place (Keating et al., 2013). In this study field scale multiphase flow reactive transport modeling was undertaken incorporating field and laboratory observations into the modeling code FEHM (finite-element-heatand-mass-transfer) (Bower and Zyvoloski, 1997). The main focus of this investigation was to explore mobilization of U, previously identified as a key risk at Chimayo. Simulations including cation

exchange/adsorption and dissolution/precipitation of calcite containing trace amounts of U were able to reproduce measured ranges and trends for pH, PCO2 , Ca, total C, U and Cl. Possibly the most integrated and applicable study to date was conducted in Mississippi (Trautz et al., 2013) and forms one of the longest and most comprehensive controlled field CO2 leakage experiments performed. In this study a CO2 charged groundwater dipole system was installed in a confined, silicate dominated aquifer over 5 months with water chemistry evolution monitored. Subsequently a reactive transport model (TOUGHREACT V2) was formed able to accurately simulate observed results by combining mineral dissolution/precipitation, cation exchange and surface complexation processes. Results showed a distinct temporal pulse signature of increased dissolved ions subsequently returning to just above background levels for most dissolved elements. This pulse effect was likely caused by rapid dissolution and depletion of reactive minerals present coupled to ion exchange/desorption. Here, a study is described which undertakes geochemical reactive transport modeling of a relatively large scale CO2 contamination field simulation conducted in spring and summer 2012 at Vrøgum plantation, western Denmark (Cahill et al., 2014). The aim of the current study is to further elucidate the risks to water resources posed by shallow aquifer CO2 contamination and determine how best to detect a leak. Furthermore, based on the work described a brief approach is suggested which can be employed to assess the risks to water resources posed by leakage and aid design of MMV programs. The current modeling study offers new insights into the impacts of CO2 contamination on a shallow aquifer system as it utilizes in-situ data obtained during controlled gas release into an aerobic, carbonate free/poor aquifer with extremely low buffering capacity, with which to calibrate the model. This type of hydro-geochemical system has not before been investigated in detail. 2. Vrøgum Investigations 2.1. Shallow aquifer CO2 contamination field laboratory The Danish research project “Environmental Technology for Geological Storage of Carbon Dioxide” (http://co2gs.geus. net/index.shtml) has been underway since 2010 exploring various aspects of CCS technology. One main focus has been the conduction of a shallow aquifer sustained CO2 contamination field experiment supported by extensive laboratory analysis. A pilot experiment was

Fig. 1. (a) Map of Denmark with field site location shown, (b) a generalized geological profile, and (c) a close up map showing the main features and the groundwater flow direction.

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

129

Table 1 Glacial sand key properties determined from field and laboratory investigations. Total inorganic carbon (TIC), quartz, K-feldspar, plagioclase and calcite are weight% composition determined by XRD. Na, Mg, Ca, K, Al and total CEC are exchange species as meq/100 g. And average background water chemistry of glacial sand used as initial model solution (mmol/l). Depth (m)

TIC (%)

Quartz (%)

Kfeldspar (%)

Plagioclase (%)

Average hydraulic conductivity (m/s)

Groundwater advective velocity (m/s)

Porosity

Na CEC

Mg CEC

Ca CEC

K CEC

Al CEC

Total CEC

06–12

0.01

94.7

1.4

3.8

2.30E-04

9.61E-07

0.3

0.01

0.15

1.08

0.03

0.04

1.31

pH 5.5

Alkalinity 0.078

O2 0.422

Cl0.943

SO4 0.16

NH4 + 0.00036

Al 0.0013

Ba 0.00017

Ca 0.091

Fe 0.0057

K 0.0597

Mg 0.168

Mn 0.00326

Si 0.129

first conducted in October 2011 followed by the main experiment, spring to autumn 2012. For full details of the site settings and experimental results the reader is referred to the pilot (Cahill and Jakobsen, 2013) and main release experiments (Cahill et al., 2014). However, for completion an overview of the field studies are given here. All field studies were conducted at Vrøgum, a managed forest plantation 3 km from the Danish west coast (Fig. 1a). Significant sediment laboratory characterization and studies were also carried out in parallel to the field investigations. The field site comprises a grassy clearing approximately 165 m by 40 m within which is a small ephemeral pond (Fig. 1c). Lithologically (Fig. 1b) the site is comprised of a 5–6 m layer of aeolian sand (including 0.2 m of topsoil) underlain by approximately 5–6 m of glacial sand with some gravel. Underlying the two upper layers, a 50 m thick formation of medium marine sand is present followed by clay at around 60 m depth. Studies have focused on the upper part (0–12 m) of the unconfined sand aquifer (i.e., aeolian and glacial sands). The water table in this phreatic formation is shallow (approximately 1.5–2 m below ground level) and the hydraulic gradient small (0.0014) with a flow direction varying between north western and south eastern area domains (from 195◦ to 185◦ , respectively) (Fig. 1c). Hydrochemistry of the Vrøgum aquifer reflects the inert silicate dominated geology; i.e., low TDS (20–30 mg/l) dominated by Na, Mg, K, Si, Ca with little or no buffering capacity (i.e., the ability to resist acidification by dissolution or precipitation of a mineral species). Distinct differences in chemistry are observed with depth and between geologic units, i.e., between the aeolian and glacial sands, likely a function of sediment/groundwater age and mineralogical composition. In the current study, modeling investigations will be in relation to the glacial sand only (reasons for which are outlined in Section 2.3). The glacial layer is a medium to coarse sand with some gravel, extremely little or no carbonate (less than 0.01% TIC) dominated by quartz, plagioclase and small amounts of K-feldspar. A summary of sediment characteristics and background average water chemistry for the glacial sand determined by standard methods during the field studies and to be employed in the model are given in Table 1. For details on methods used for determination the reader is referred to the pilot and main injection experiments previously cited.

2.2. Batch reactor experiments During field site characterization and experimental design, a series of batch reactor experiments were conducted to elucidate likely aqueous chemistry changes induced by injection of CO2 in the field. During the experiment 50 g aliquots of glacial sand (obtained from coring to 10.2 m depth) were weighed in duplicate into 300 ml glass flasks with 100 ml of de-ionized water, capped and sealed with a rubber septum and equilibrated for seven days. Similar batch reactor experiments (unrelated to the current field site) observed relatively stable background chemistry after 6 days pre-CO2 incubation for a range of mixed composition sediments (Lu et al., 2010;

Na 0.857

Cahill et al., 2013). After the equilibration period background samples (10 ml removed by syringe) were taken with pH determined immediately (Hach Lange probe, accuracy 0.2) and dissolved ion concentrations by ICP–MS analysis (PerkinElmer Elan 6100 DRC, accuracy 5%) following preservation (filtration with 0.45 ␮m cellulose acetate filters and acidification to pH 2 with concentrated HNO3 ). Following the initial background sample, batch reactors were exposed to food grade CO2 for 30 min via a needle injection and venting system. CO2 was injected directly into the liquid and sediment mixture causing turbulent mixing at an overpressure of approximately 1.5 bars venting through the head space. After CO2 gas injection, overpressure was allowed to equalize with atmospheric pressure and the reactors were closed. Samples were subsequently taken after 13 days with pH and dissolved ions measured as described for the background samples. This methodology allows the impacts of CO2 saturation on reactor water chemistry to be determined and is not intended to perfectly mimic or predict field conditions. The premise for these reactor experiments is that a pre- and post- CO2 batch reactor water chemistry is determined in a highly controlled, simple and stagnant system. This allows focus on the CO2 induced geochemical mechanisms alone (i.e., complications of flow and heterogeneity are minimized) and indicates what will likely happen in the field.

2.3. Field study results Following site characterization, laboratory sediment studies (including the batch reactor studies previously outlined) the main field release experiment was designed and constructed. During the main release experiment 1600 kg of gas phase CO2 was delivered to the subsurface through 4 inclined injection wells (2 in the aeolian and 2 in the glacial sand) over 72 days. The resultant dissolved plume of CO2 was subsequently monitored over 252 days and 20 m of flow with high temporal and spatial resolution. Physico chemistry (pH, EC and dissolved oxygen by Hach Lange HQ40D meter, intelliCALTM probes and flow cell), alkalinity (Gran titration method) and all aqueous cation concentrations (ICP–MS: PerkinElmer Elan 6100 DRC) were determined during every sample round. Anion concentrations were determined for only several sample rounds. This method was justified based on results seen in other studies where it was seen that anion concentrations (other than HCO3 − ) show little or no change as a result of CO2 dissolution (Peter et al., 2012; Trautz et al., 2013). PHREEQC calculations of saturation indices of groundwater chemistry using measured compared to average anion concentrations, showed very small differences, confirming this is a reasonable simplification. Therefore for ease of calculation input, a fixed average background concentration method was employed. For more information see Cahill et al., (2013). Results for the main experiment showed evolution of 2 seemingly distinct plumes in the aeolian and glacial sands. Dissolved CO2 plumes showed complex and distinct evolution of physico and aqueous ion chemistry with flow and distinct behavior of dissolved ions observed and categorized into groups. The glacial

130

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

Fig. 2. (a) Center line cross section of Vrøgum groundwater saturation for CO2(gas) (SI = log(PCO2 ) after 52 days of injection, (b) cross section of electrical conductivity within high density sample line normal to flow 1.5 m down gradient of injection plane after 29 days of injection, and (c) electrical conductivity at 8 m depth in plan view at various times pre and post injection commencement. Sample points denoted by +(measured during sample round) or ♦ (not measured during sample round).

sand layer (approximately 6–12 m depth) was seen to be relatively homogenous (i.e., in terms of mineralogy and physical properties) and showed the most consistent plume development in comparison with the more heterogeneous aeolian sand. A uniform plume of dissolved CO2 was observed to evolve with flow and move through the formation in a very stable manner moving consistently through the sample points (Fig. 2). For this reason the glacial sand was chosen as the modeling target in order to allow for a simplified 1D model. 3. Geochemical model development 3.1. Conceptual model Due to the injection system setup and the relatively high homogeneity of the glacial sand (in terms of mineralogy and hydraulic properties) a very simple conceptual model of the system can be formed on which a model can be based. The glacial sand is essentially a column experiment within a natural system. CO2 is introduced at the start of the column forming a reasonably homogeneous vertical plane of elevated PCO2 (Fig. 2). CO2 charged groundwater then moves with flow through multiple sample points providing data along the full length of the column. Given the width of the plume and that in groundwater flow systems such as this, transverse dispersion is negligible (Freyberg 1986; Jensen et al., 1993; Engesgaard et al., 1996) it is reasonable to model the system in 1D. PHREEQC (Parkhurst and Appelo, 2013) is a geochemical modeling program from which simple 1D transport models have been successfully developed previously (Postma and Appelo, 2000; Postma et al., 2007). Reactive transport simulations can include the full range of processes likely occurring in the glacial sand and comprise a series of cells with specified sediment properties and mineral phases arranged in a column through which flow occurs starting with an initial solution. Each transport step involves the replacement of initial solution with the defined solution in cell 0 and implementation of the geochemical processes

desired. Simultaneously all other cell solutions are shifted “down flow” and enter the geochemical processes imposed in the given cell. Hydrodynamic dispersion is handled by mixing adjacent cells. In PHREEQC the field system is simulated by cells containing initial solutions of average glacial sand water chemistry (as 1 L of pore water) with an inflow solution (also equivalent to the background chemistry) charged with CO2 gas (by imposing equilibrium with gas phase CO2 at an applicable partial pressure) being continuously introduced. Processes are then included within the cells to simulate the chemistry changes observed. A schematic of the PHREEQC 1D reactive transport model concept is given in Fig. 3. Average glacial sand pore water chemistry used as the inflow and initial cell solution determined during background monitoring is given in Table 1. 3.2. Batch and reactive transport models Based on laboratory analysis, field results and other published studies, e.g., (Zheng et al., 2012; Trautz et al., 2013), geochemical processes which may be controlling CO2 induced water chemistry change in the glacial sand aquifer were identified and are listed in Table 2. Ion exchange was included based on laboratory determined cation exchange capacity of 1.3 meq/100 g sediment and considered in the model, both by using the standard phreeqc.dat ion selectivity coefficients and using fitted coefficients. Fitted selectivity coefficients were based on the observed natural water chemistry composition and the in-situ exchangeable cation composition – fitted using PHREEQC in batch mode. A goethite surface and gibbsite presence/abundance were included based on Fe and Al release for sediment extraction results using an NH4 -oxalate and ascorbic acid extractant (which targets crystalline, hydrous Fe and Al oxides). A ferrihydrite surface and amorphous gibbsite presence/abundance were included based on Fe and Al release for sediment extraction results using an NH4 -oxalate extractant which targets amorphous and poorly crystalline hydrous Fe and Al oxides. Inclusion of pH dependent dissolution of amphibole and feldspar (i.e., lower pH

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

131

Fig. 3. Model concept employed in PHREEQC. Series of cells hosting geochemical processes through which average background chemistry of glacial sand groundwater is shifted. Red dots represent potential points where direct field experiment measurements are available for comparison. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 2 Geochemical processes investigated during model formation. Geochemical process

References

Aqueous complexation Ion exchange (standard exchange coefficients) Ion exchange (fitted exchange coefficients)

wateq4f.dat PHREEQC database phreeqc.dat PHREEQC database Fitted in batch mode to measured exchanger composition based on background water chemistry (Jessen et al., 2012) (Dzombak and Morel, 1987) in phreeqc.dat phreeqc.dat phreeqc.dat Formed for this study Formed for this study phreeqc.dat

Goethite surface complexation Ferrihydrite surface complexation Gibbsite equilibrium control Kaolinite equilibrium control Kinetically controlled dissolution of gibbsite pH dependent amphibole/feldspar dissolution rates Ferrihydrite equilibrium control

Final model composition √ × √ × × √ × × × ×

Table 3 Measured exchanger composition, PHREEQC.dat and fitted exchange coefficient exchanger composition and PHREEQC.dat and fitted exchange coefficients (both shown relative to Na+ following the Gaines-Thomas convention). Exchange species

Average measured glacial sand exchanger composition (meq/l pore water)

Calculated exchanger composition using phreeqc.dat exchange coefficients (meq/l pore water)

Exchanger composition from fitted selectivity coefficient (meq/l pore water)

phreeqc.dat selectivity coefficient

Fitted selectivity coefficient

NH4 + Na+ Mg2+ Ca2+ K+ Al3+ Ba2+ Sr2+ Zn2+

0.22 2.18 9.71 69.20 1.74 2.33 0.10 0.05 0.05

0.01 2.00 44.28 40.06 0.92 1.99 0.09 0.28 0.16

0.27 3.59 8.41 72.62 2.08 2.25 0.07 0.03 0.04

0.25 1.00 0.50 0.40 0.20 0.73 0.35 0.35 0.40

0.03 1.00 2.60 0.60 0.16 1.08 0.75 1.78 1.50

incurs more dissolution) was considered during initial model set up (as Mg and Ca were key ions of interest in the system). Moles of silicate phases were added (i.e., dissolved) stepwise as a function of pH, if pH was below 5.1 (based on observation made by (Sverdrup and Warfvinge, 1988)). The goethite and ferrihydrite surfaces were defined according to (Jessen et al., 2012) and (Dzombak and Morel, 1987) respectively. Cation selectivity coefficients (phreeqc.dat and fitted) including subsequent modeled batch mode exchanger compositions are given in Table 3. Surface complexation parameters used in the model and abundance of hydrous Al and Fe oxides present in the sediment are given in Table 4.

Batch reactor results were simulated in PHREEQC in batch mode by simulating experimental conditions assuming a number of equilibrium processes. A 100 ml solution of measured background sample composition was used as input with combinations of geochemical processes included, however exchange and surface complexation parameters were modified to represent the different sediment water ratio in batch reactors (i.e., only 50 g of sediment present). Numerical batches were now subjected to imposition of 1 bar CO2 (as CO2 gas saturation index of 0 in PHREEQC) with inclusion of gibbsite (amorphous and crystalline) equilibrium, cation exchange and surface complexation in different combinations to

Table 4 Surface complexation parameters and mineral quantities. Amorphous and crystalline Al and Fe oxides are moles/kg sediment determined by targeted sediment extractions. Surface Goethite Ferrihydrite Mineral Quantities Amorphous hydrous Al oxides 0.002

Surface area (m2 /g)

Amount (g/kg)

Uni: 3.4 Tri: 2.8 (per nm ) Weak: 0.023 strong: 0.0005 (moles)

60 600

5.8 3

Crystaline hydrous Al oxides 0.007

Amorphous hydrous Fe oxides 0.004

Crystaline hydrous Fe oxides 0.066

Sites 2

132

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

1.6E-04

1.8E-04 1.6E-04

1.2E-04

Concentraon (Moles/l)

Concentraon (Moles/l)

Goethite Surface

2.0E-04

Exchange

1.4E-04

1.0E-04 8.0E-05 6.0E-05 4.0E-05

1.4E-04 1.2E-04 1.0E-04 8.0E-05 6.0E-05 4.0E-05

2.0E-05

2.0E-05

0.0E+00

0.0E+00

1.8E-04

Ferrihydrite Surface

1.6E-04

Exchange and Gibbsite

1.6E-04 Concentraon (Moles/l)

Concentraon (Moles/l)

1.4E-04 1.2E-04 1.0E-04 8.0E-05 6.0E-05 4.0E-05

1.4E-04

Pre CO2

1.2E-04

Post CO2

1.0E-04

Modelled

8.0E-05 6.0E-05 4.0E-05

2.0E-05

2.0E-05

0.0E+00

0.0E+00

Ca

Exchange and Goethite Surface

2.0E-04

Mg

Na

K

Al

Si

Fe

Mn

Sr

Ba

Zn

6.0

1.8E-04

1.4E-04 5.0

1.2E-04 pH

Concentraon (Moles/l)

Batch Reactor: Actual and Modelled pH

5.5

1.6E-04

1.0E-04

4.5

8.0E-05 6.0E-05

4.0

4.0E-05 2.0E-05

3.5

0.0E+00

Ca

Mg

Na

K

Al

Si

Fe

Mn

Sr

Ba

Zn

Pre CO2

Post CO2

Exchange

Exchange Goethite Ferrihydrite Exchange and Gibbsite Surface Surface and Goethite

Fig. 4. Batch experiment pre and post CO2 chemistry and batch model results showing output for varying geochemical process inclusion.

determine if a good approximation to results could be attained. Batch model results are shown in Fig. 4. A basic 1D flow model was formed in PHREEQC comprised of 50 cells and transport steps, flux flux boundary conditions (at each end of the forward flowing model) and dispersivity and diffusions coefficients of 0.1 m and 0.5E–9 m2 /s, respectively. Each time step moves water from 1 cell to the next in 2.5 days and as such cell lengths (inferring advective velocity) were calibrated to fit breakthrough of dissolved CO2 as detected by EC, pH and dissolved ion changes during the field experimental phase. It was seen that advective velocity was slower (approximately 0.8 m/day) in the first few meters from the injection plane increasing (to 0.16 m/day) approximately 1.5 to 2 m down flow. Again geochemical processes which may be controlling CO2 induced water chemistry change in the glacial sand aquifer (listed in Table 2) were considered in various combinations. This methodology allowed the assessment of potential contribution to water chemistry change for each process and the most applicable model in terms of being as simple as possible, but still able to describe the observations adequately to be formed. Field results and a selection of the best obtainable model outputs of varying process inclusion for 35 days post injection commencement with distance from the injection plane along the center line for pH and major ions is given in Fig. 5. Results for pH dependent dissolution of silicates (amphibole and feldspar) are not shown however are briefly discussed in the proceeding section. Introduc-

tion of CO2 into the flow model (i.e., simulation of gas injection) is achieved by imposing a saturation of the gas (i.e., partial pressure) on the inflowing solution. 3.3. Geochemical process contribution and selected model From Figs. 4 and 5 individual contributions of various geochemical processes can be seen. From these simulations several key points regarding geochemical model construction can be construed; (1) cation exchange (using fitted selectivity coefficients) and a Goethite surface contribute to an increase in dissolved ions alone inferring dissolved CO2 can induce desorption from the exchanger and mineral surfaces. In the model this is due to aqueous complexation whereby CO2 charged conditions changes element activities in solution, e.g., more Al is present as Al3+ , thus small changes in exchange and surface bound ion exchange are induced. Alone, these surface based processes vastly exaggerate pH decreases. (2) Exchange and gibbsite equilibrium vastly improves pH and major/minor ion simulations. (3) Inclusion of a goethite surface with exchange provides little to no additional benefit to the model fit. Presumably specific effects of surface complexation related to pH changes were small enough that it was not necessary to include in the model. Additionally, computational time is increased dramatically by inclusion of this process with some parameters inaccurately fitted at later times. 4. Inclusion of a limited pool of

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

133

Fig. 5. Reactive transport model output for varying geochemical process inclusion in comparison to field results with distance after 35 days of CO2 injection along the centerline in the glacial sand.

reactive, amorphous gibbsite enables the initial higher concentrations of all ions at the plume front to be simulated. The inclusion of two equilibrium controlled gibbsite phases is made possible by the dissolution only option for equilibrium phases in PHREEQC. As intuitively expected, inclusion of pH dependent dissolution of silicate minerals, without kaolinite equilibrium imposed was seen to give Si concentrations proportional to the level of silicate added. However it did not improve the fit for the cations, which was the intention. Furthermore the most important observations could be simulated in the chosen model with fewer processes and fit parameters, by parameters for which clear evidence for their influence was attained (e.g., measured parameters such as CEC composition). Consequently silicate dissolution was not included in the further development of the model and results are not presented. Limitations associated with this decision are discussed further in Section 3.5. From trial simulations, imposition of a CO2(gas) saturation index (SI) of −0.15 (corresponding to a PCO2 of 0.7 atm, within the range observed in the field) within a 1 dimensional flow model including; 1. aqueous complexation, 2. ion exchange with fitted selectivity coefficients and 3. Limited (i.e., finite quantities) amorphous and crystalline gibbsite equilibrium (dissolution only) – was able to simulate most observed field measurements reliably with distance and time (Figs. 6 and 7). In summary the gibbsite and amorphous Al-hydroxide phases are used as the reactive mineral phases that dissolve (i.e., during buffering) to maintain saturation equilibrium, while the ion exchange in the model is assumed to take place on an unspecified Al–Si-layer-silicate, such as a montmorillonite, with a fixed total exchange capacity independent of the gibbsite. Major (Ca, Mg, Na, K) and minor (Al, Ba, Sr) ions, pH, Alkalinity and PCO2 (as SI) were all approximated reliably over the greatest distance and time scales with this model composition. Field and final model results with time (i.e., breakthrough) at 1, 2 and 5.25 m down flow and results with distance for 35, 65 and 79 days are shown in figures 6 and 7, respectively. With groundwater observed moving approximately 0.08 m/day through the injection plane, it is reasonable to assume that displacement of CO2 charged groundwater has not yet begun in the first down-gradient sample points (i.e., >1 m down flow) after 7 days. Consequently despite cessation of injection after 72 days it is assumed the effects of constant CO2 injection still dominate after 7 days, justifying inclusion of the 79 day sample. Some parameters were not successfully simulated using any combination of processes in the model, i.e., Zn, Mn and Fe the potential reasons for which will be discussed in the proceeding sections. Prior to CO2 injection simulations, a model simulation was run under ambient

conditions (i.e., with average natural PCO2 levels throughout) and was seen to provide stable background chemistry throughout the model domain. This stable base model will allow the impacts of CO2 on water chemistry to be modeled and investigated further. 3.4. Model sensitivity and mechanisms Sensitivity of the final flow model to various parameters was assessed during development with parameters identified to elicit significant control on modeled results. PCO2 (i.e., CO2 saturation index) employed, as would be expected, exerts significant control on the severity of the response in terms of pH decreases and modeled ion concentration increases. This parameter however is relatively well constrained, showing only small fluctuations along the flow line during the experimental period, so, an overall SI of −0.15 (0.7 atm) as used in the model for the in-flowing solution is justified. Saturation indices imposed for equilibrium control by amorphous and crystalline gibbsite were also seen to exert some control on results. By increasing or decreasing the saturation index (corresponding to a change of the log K value by the same amount) imposed above or below 0 greater or lesser responses in terms of pH decrease and dissolved ion increases can be simulated. Overall, setting an SI at or around 0 for these minerals allows the most accurate simulation of modeled parameters inferring control of water chemistry evolution by gibbsite equilibrium is a key process at Vrøgum. Control of groundwater chemistry by gibbsite equilibrium including an amorphous form has been previously inferred in aquifer acidification studies supporting this hypothesis (Dahmke et al., 1986; Reuss et al., 1987; Hansen and Postma, 1995; Kjøller et al., 2000). Furthermore these studies also observed small changes in aqueous Al related to much larger changes in Mg and Ca concentrations as seen in the current study. For more discussion on and support for this proposed mechanism the reader is referred to the full field experiment article (Cahill et al., 2014). The amount of gibbsite like minerals included in the model also determined model fit significantly. A limited amount of the amorphous form, i.e., 0.00016 moles per cell (i.e., 0.002 g/kg of sediment), was seen to create a plume front peak of amplitude and length in concurrence with observed results. This amount is rather low in comparison to the amount potentially present determined by NH4 -oxalate targeted extractions of 0.11 g/kg. This likely reflects differences in model operation, in-situ geochemical processes and the extraction method. In the model all mineral present is permitted to react to equilibrium (i.e., full aqueous contact is assumed). In reality under

134

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

Fig. 6. Modeled and measured (a) major and (b) minor ion concentrations vs. time with distance from the injection plane generated from final selected model including cation exchange with coefficient fitted to batch measurements combined with gibbsite (amorphous and crystalline) equilibrium.

field conditions not all reactive amorphous gibbsite may be available to react (i.e., in aqueous contact) and the pool of mineral will not be simply split between amorphous and crystalline with a continuum of stabilities more likely, as demonstrated by (Larsen et al., 2006) for Fe-oxides. Furthermore, the extraction utilizes agitation/mixing and a pH of 3 (compared to around 4.5 observed and modeled) to ensure complete reaction and also targets Fe oxides in which Al may be a substitute cation. Thus the extraction method should indicate a larger amount of gibbsite is present than required in the model and it is reasonable to assume less amorphous gibbsite

a

Ca

pH measured

Ca measured

Mg

Na

K

Mg measured

Na measured

K measured

65 Days

35 Days

2E-03

5.5

79 days

5 2E-03

4.5

pH

Concentraon (moles/l)

pH

would be required to simulate field results. The amount of crystalline gibbsite included in the model determines the secondary peak length for major and minor ions and the extent of acidification, the implications of which will be discussed in the proceeding sections. Cation selectivity coefficients exert some influence on model output as described in Sections 3.3. and 3.4. The fitted selectivity coefficients (i.e., fitted against the natural exchanger composition and pore water concentrations) do give a better fit compared to the standard phreeqc.dat coefficients. The allowance of dissolution and precipitation of gibbsite also exerts an effect on modeled results.

1E-03

4 5E-04

3.5

3

0E+00 0

1

2

3

4 Distance (m)

5

6

8 0

7

b

Al

5

Sr

10 Distance (m)

Ba

Al measured

Concentraon (moles/l)

1.2E-05

15

20 0

Sr measured

5

10 Distance (m)

15

20

Ba measured

65 Days 79 days

35 Days

1.0E-05 8.0E-06 6.0E-06 4.0E-06 2.0E-06 0.0E+00 0

1

2

3

4 Distance (m)

5

6

7

8

0

2

4

6

8

10 12 Distance (m)

14

16

18

20

0

2

4

6

8

10 12 Distance (m)

14

16

18

20

Fig. 7. Modeled and measured (a) major and (b) minor ion concentrations vs. distance along the centerline with time generated from final selected model including cation exchange with coefficient fitted to batch measurements combined with gibbsite (amorphous and crystalline) equilibrium.

18

135

57

CaX

16

MgX

NaX

KX

AlX

SI CO2

14

56

12 55

10 8 6

54

4 2

Ca composion (meq)

exchanger composion (meq) and CO2 saturaon index

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

53

0 52

-2 0

2

4

6

8 Distance (m)

10

12

14

16

Fig. 8. Model exchanger composition (meq/l pore water) with distance along the centerline including modeled saturation index for CO2(gas) after 65 days injection.

When both processes are included, a sufficient buildup of Al is not achieved because the amorphous gibbsite is immediately turned into the more crystalline gibbsite, which in turn causes underestimation in displaced cations compared to observed field results. Imposition of “dissolution only” allows field results to be simulated accurately using equilibrium control of two types of gibbsite. This requirement does not necessarily imply that precipitation is not occurring in the field; rather the rate of dissolution is faster than precipitation under the conditions imposed by CO2 . Total cation exchange capacity appears to have a small to moderate effect on model results, i.e., changes of 20% have little effect on observed results. Based on conducted simulations, mechanisms for simulating CO2 induced water chemistry changes in the final model can be further described by considering the composition of the exchanger as shown in Fig. 8. The model simulates the observed field results by gibbsite-derived Al3+ driven exchange processes, originating from both an amorphous (i.e., more reactive) and crystalline (i.e., less reactive) pool of the mineral. Reaction equations and equilibrium constants for these Al-hydroxides are given in Table 5. As a result of this gibbsite based equilibrium control, the modeled exchanger surface shows an increase in exchanged Al starting in the first cell, which becomes depleted of other elements in turn, based on affinity for exchanger and abundance. This process proceeds cell by cell the extent of which is dependent on the amount of crystalline gibbsite included. As a cells exchanger takes up Al and other elements are displaced the adjacent down flow cell initially shows a small increase for displaced elements (particularly Ca) reflecting the increased concentrations in solution and abundance on the exchanger. Assuming an infinite amount of crystalline gibbsite, cell exchangers in the model would eventually all reach a steadystate controlled by the background concentration of other cations entering from upstream, the PCO2 , controlling the pH, which again determines the Al concentration determined by equilibrium with gibbsite.

3.5. Model limitations Although considering many geochemical processes in various combinations, the final flow model was unable to simulate observed results for Mn, Fe, Zn and Si. Mn and Fe were both observed to decrease following injection of CO2 into the glacial sand creating an apparent inverse plume of imposed zero concentration with flow. This effect was potentially attributed to precipitation into a Fe-bearing Al-hydroxide, however no conclusive evidence of this was found, more discussion upon which can be found in the field experiment description paper (Cahill et al., 2014). Another possibility is that colloidal fractions play a role in evolution of these elements aqueous chemistry’s however further investigation of this is beyond the scope of the current study. Zn was generally seen to be far in excess of modeled predictions based on Al induced ion exchange alone inferring another source of Zn within the sediment. A Zn bearing mineral was not detected during any of the sediment analysis conducted inferring maybe trace amounts (i.e., below levels of detection) of a dissolving Zn-bearing mineral were present or Zn was present as a substituting cation in another dissolving mineral, neither of which hypotheses can be proven with the available data. Slight discrepancies were seen at certain spatial and temporal points between the model and field data in particular for trace metals and pH, the cause of which cannot be determined. However, localized heterogeneities (e.g., local pockets of increased gibbsite/calcite and/or silt) are potentially present within what is otherwise relatively homogenous sand. As mentioned, including Si weathering did not improve the fit for cations, implying that fitting the Si was pure fitting. Si concentrations could however be simulated reasonably in the flow model by inclusion of a goethite surface however this process made the model inefficient (i.e., significantly increasing run time) and unstable (i.e., causing convergence problems). As Si is not an element of major concern with regards to water quality deterioration and was present only in low concentrations it was omitted from the final model. This omission however remains a

Table 5 Key equations involved in model operation including; Al-hydroxide dissolution (a gibbsite and amorphous type), and Al3+ -driven ion exchange (shown for Al–Ca exchange) and associated standard log K values and imposed saturation indices. Please note that imposing a change in the saturation index effectively changes log K by an equivalent amount. Reactions Gibbsite dissolution Amorphous Al–OH dissolution Al3+ exchange

Equations +

3

Al(OH)3 + 3H ↔ Al + 3H2 O 3 Al(OH)3 + 3H+ ↔ Al + 3H2 O 3+ 1 1 1 Al + ca − X ↔ Al − X3 + 2 3 2 3

1 ca2+ 2

Log K

Imposed saturation index

8.1 10.8 -0.8

-0.1 0.2 n/a

136

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

Fig. 9. Future predictions using calibrated reactive transport model for (a) major ions concentrations and pH and (b) Al concentration. The 3 year major ion and pH graph shows the generalized phases of leakage as identified in the field and modeling study: advanced acidification (A), post-pulse phase (B) and pulse phase (C).

limitation as changes in concentrations of Si are relevant to silicate dissolution and secondary mineral formation and transformations which could be important processes in a system such as Vrøgum, especially at much longer timescales. Following development of a model able to simulate the observed field results to an acceptable degree, a series of model predictions were made including several scenarios. It should be kept in mind that the model is not calibrated for these scenarios, so model results should be seen as probable outcomes. Modeling related to carbon sequestration is a field with many limitations and challenges, described in several review papers on the topic, e.g., (Gaus et al., 2008).

4. Predictions and scenarios 4.1. Forward projection The selected model was extended spatially and temporally in order to predict the likely effects on water chemistry and sediment at Vrøgum over timescales applicable to CCS leakage (i.e., several years or longer). Fig. 9 shows water chemistry for modeled parameters after 1–3 years. Predictions show clearly the pulse of increased ions moves with advective flow behind which aqueous element concentrations falls in two steps with zones of increasing acidification formed. Key characteristics are the zone of maximum acidification (caused by buffering mineral exhaustion) behind the initial pulse (induced by amorphous gibbsite dissolution) and stable secondary phase (induced by crystalline gibbsite dissolution) which grows with time at a rate slower than advective flow. Growth of the zone of advanced acidification is controlled by the time needed to deplete all forms of gibbsite (the sole and primary buffering mineral species in this specific case). Al concentration evolution follows other dissolved elements with an initial pulse followed by an elevated post-pulse phase. However, a pronounced secondary peak is now present before decrease to background levels. The secondary Al spike occurs at the interface of buffering exhaustion (i.e., where fully depleted cells meet gibbsite containing cells) and

is caused by dispersion of inflowing CO2 charged, gibbsite subsaturated acidified groundwater. 4.2. Depth increase The model was altered to reflect greater depth by increasing SI of CO2 (inferring higher PCO2 ) applied to the inflow solution (representing 2, 5 and 10 bar) in order to assess the potential effects over a range of depths which may be impacted by leakage of CO2 . Results for 65 days of simulation are shown in Fig. 10. Results infer increasing depth induces greater acidification upon buffering mineral exhaustion, i.e., pH of 3.8, 3.6 and 3.4 for 2, 5 and 10 bar PCO2 , respectively. Also the rate of aquifer acidification is more rapid with depth (due to faster depletion of the buffering mineral, i.e., gibbsite) with the zone of advanced acidification seen to extend approximately 3.5, 6.5 and 7.5 m from the leakage plane after 65 days simulation for 2, 5 and 10 bar, respectively. Concurrently Al concentrations are seen to increase with depth (i.e., increasing PCO2 ) exhibiting greater but less extensive peaks. The secondary peak of Al at the buffering exhaustion interface is also still evident as described in Section 4.1. 4.3. Calcite inclusion Calcite was separately introduced to the model (using observed average field SI, i.e., PCO2 fixed at −0.15) to assess the effects of varying buffering capacity on water chemistry evolution (0.1 and 0.5% calcite sediment composition by dry weight) as an equilibrium phase. Results for 65 days of simulation are shown in Fig. 11. Calcite equilibrium has previously been seen to dominate water chemistry in CO2 contaminated water systems (Cahill et al., 2013). Adding a new reactive phase (i.e., calcite) which was not included in the construction of the model should be viewed only as a demonstration of the potential effects of the presence of carbonates (i.e., buffering) rather than a realistic simulation. Addition of calcite reduces the level of acidification observed significantly. No difference in water chemistry evolution for calcite levels greater than 0.5% was apparent (hence only 0.5% is shown) inferring even small amounts of the

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

137

Fig. 10. Effects of increasing PCO2 on (a) major ions concentration and pH development and (b) Al concentration development using constructed reactive transport model after 65 days of simulated CO2 injection.

a

Alkalinity

Ca

Mg

Na

K

b

pH 1E-04

7.5

0.5 % Calcite

7

2E-02

6.5 6

2E-02

5.5 5 1E-02

4.5 4

5E-03

concentraon (moles/l)

0.1 % Calcite

pH

concentraon (moles/l)

3E-02

Calcite Content 0.1% 0.5%

8E-05

6E-05

4E-05

2E-05

3.5 3

0E+00 0

5

10

15

distance from injecon (m)

20

0

5

10

15

distance from injecon (m)

20

0E+00 0

5

10

15

20

distance (m)

Fig. 11. Effects of varying calcite content (with a fixed PCO2 inferred as saturation index of −0.15) on a) major ions concentration and pH development and (b) Al concentration development using validated reactive transport model after 65 days of simulated CO2 injection.

mineral provide significant buffering capacity. Although limiting acidification, calcite severely increases Ca concentrations and also other aqueous ions (via increased Ca-driven ion exchange). Sediment with only 0.1% calcite composition initially displays the same level of buffering as more carbonaceous formations however this capacity is quickly exhausted. Upon exhaustion of calcite in the first few cells of the model domain, pH in this scenario falls further. Al concentrations for model scenarios with greater than 0.5% calcite remain low and relatively unchanged as the high buffering capacity prevents significant mobilization. For the 0.1% calcite scenario, once depleted pH falls further, Al is seen to increase rapidly as gibbsite now dissolves more readily in this zone of advanced acidification. 5. Discussion 5.1. Geochemical process contribution During model development, the contributions of individual geochemical processes to water chemistry change caused by CO2 introduction were explored. A combination of gibbsite equilibrium (dissolution only) coupled to ion exchange (with Vrøgum specific fitted selectivity coefficients) was required to simulate a response similar in magnitude to that elicited by CO2 in the glacial sand.

The exact phase responsible for the exchange processes was not identified in the current study however based on XRD results a small clay mineral fraction most likely provides the exchange surface. In the selected model selectivity coefficients were fitted so that the exchanger composition simulated by PHREEQC from measured pore water concentrations better matched the sediment’s measured exchanger composition (Table 4). This method is not a calibration of selectivity coefficients to fit simulated results to observed experimental field concentrations, but calibration of the exchange system prior to simulations to represent the Vrøgum glacial sand. The fitted selectivity coefficients do deviate from published values, however not by an alarming degree (i.e., within or close to ranges reported in the literature) as selectivity coefficients have been reported over considerable ranges (Bruggenwert and Kamphorst, 1979). Thus it is possible that Vrøgum’s glacial sand is different from other investigated sediments and variation from other published results should not infer an issue. A comparison of values used in this study and those employed by the PHREEQC.dat database (i.e., published values) is made in Table 3. For example PHREEQC.dat uses values determined by (Laudelou et al., 1968) for Ca and Ba with selectivity coefficients of 0.4 and 0.35. The current study observed values of 0.6 and 0.75 for Ca and Ba respectively were more appropriate.

138

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

Inclusion of finite amounts of amorphous and crystalline gibbsite in conjunction with equilibrium maintenance and a dissolution only condition was also required to simulate most accurately observed field results. The hypothesis that forms of gibbsite of varying stability are present and under equilibrium control is logical and reasonable based on the data presented in the field study description (Cahill et al., 2014) and other cited aquifer acidification studies (Hansen and Postma, 1995; Kjoller et al., 2004). A more complex model involving a kinetic description of the dissolution of a continuum of gibbsite reactivity’s could presumably be made to fit the observations also. However, given that we know very little on the dissolution and precipitation kinetics of the gibbsite type mineral present, this would add uncertain parameters to the model and so has not been undertaken for the present study. Additionally, kinetic inclusion of other minerals present in the glacial sand which could theoretically also be involved in dissolution reactions which may influence chemistry evolution (e.g., plagioclase and amphibole) would also be desirable. However, as these minerals are significantly less reactive than gibbsite and the current model is able to explain the majority of chemistry change observed with relative accuracy within the given time, it was not explored in the current study. Though quartz, makes up approximately 95% of the glacial sand, it is well known to be highly inert and very unreactive, particularly under observed and induced conditions (i.e., low pressure, temperature and pH > 4.8) and time scales (72 days of injection). Consequently, in order to simplify model operation, quartz was not included. Overall the model formed demonstrates that Al3+ -driven ion exchange, originating from dissolution of gibbsite of varying reactivity is able to reproduce the observed chemistry reasonably well. However, this is only one realization and other model set-ups may be able to attain comparable results such is the nature of modeling. The final model should not be inferred as a definitive account of the mechanism of chemistry change at Vrøgum, rather a representation of the most important geochemical mechanisms most likely occurring based on data attained from field, laboratory and previously mentioned groundwater acidification studies (Hansen and Postma, 1995; Hansen and Postma, 1995). The discrepancy between fitted and commonly used selectivity coefficients, the use of two ‘separate’ gibbsite pools (crystalline and amorphous) and inclusion of a dissolution only condition are indications that the model is a simplification. It is possible other undetected or unidentified trace minerals are also contributing to, or causing, observed chemistry evolution. For example trace amounts of calcite may be contributing, as suggested in other studies (Zheng et al., 2012). However, calcite inclusion simulations (Section 4.3) show even small amounts of calcite would have induced different chemistry than observed (i.e., higher alkalinity, smaller pH decreases, higher Ca coupled to lower Al concentrations). In addition the flow systems conceptual and numerical model is further supported independently by batch model results whereby the same model parameters most accurately reproduce the observed sediment-water CO2 exposure results. The inclusion of only finite amounts of gibbsite of varying reactivity is logical and was seen to allow continued acidification to be simulated as observed in the field. This effect infers a sediment depletion (or buffering mineral exhaustion) effect would occur during prolonged CO2 leakage. This phenomenon has significant implications for both risks posed to water resources by leakage and for MMV program design and implementation. 5.2. Spatial and temporal chemistry evolution Field and model results clearly show CO2 induced chemistry evolution is transient and spatial, thus risks to water resources would change with time and distance from a carbon sequestration leakage horizon. Both field result simulations and predictive

model results (Figs. 6, 7 and 9) show the development of an initial pulse of increased ion concentration as seen in other studies (Cahill and Jakobsen, 2013; Trautz et al., 2013). The magnitude and duration of the pulse was seen to be directly influenced by the reactivity and amount of buffering minerals present (i.e., the minerals providing the buffering capacity; amorphous and crystalline gibbsite), the depth of leakage (input as increasing PCO2 ) and distance from the leakage feature. During this pulse phase, buffering by amorphous gibbsite takes place. Consequently Al concentration increases sharply, inducing exchange processes coupled to a small pH reduction of approximately −0.2 units. The elevated pulse phase is seen to move with advective flow becoming broader with distance as more amorphous gibbsite is encountered and dissolved. Once the pool of reactive gibbsite is depleted elemental peaks decrease (including Al), reflecting a lower rate of exchange due to less Al being released into solution. Behind the pulse phase the system is now buffered by crystalline gibbsite whereby Al concentrations slowly and steadily rise. This second phase leads to formation of a zone of aquifer characterized by further reduction in pH (an approximately −0.3 unit further decrease) and elevated but slightly lower concentrations of exchanger derived ions. Model simulations predict that eventually crystalline gibbsite is also exhausted leading to a zone of further acidification (additional reduction in pH of −0.8 units). At this point there are no mineral phases remaining in the model to buffer (i.e., be consumed whilst resisting acidification) the CO2 acidified groundwater. A buffering mineral is defined as any mineral species reactive and abundant enough that it can, by dissolution (i.e., its own consumption and that of protons), resist acidification. Following complete exhaustion of buffering minerals, exchanger derived aqueous ion concentrations return toward background levels with only a small elevation above natural conditions is still evident. This is a result of aqueous complexation as previously described. The zone of extensive acidification moves through the model slower than advective flow due to the amount of crystalline gibbsite in the sediment (i.e., a greater amount than the dissolution capacity of the incoming water). Forward projections suggest the zone of extensive acidification would expand 50 m from the leakage horizon in 3 years, presumably growing indefinitely whilst leakage continued. This is an unacceptable consequence of leakage considering many countries such as the U.S and Denmark demand no measurable impact. In the model, remaining sediment is assumed inert; however in reality this would not strictly be true. All minerals (e.g., silicates and even quartz) will buffer over appropriate timescales, theoretically providing an infinite buffering capacity. At this point less reactive silicates and other inert minerals would now dominate the observed geochemical processes occurring according to the Goldich dissolution series (Goldich, 1938). The concept of buffering exhaustion is nevertheless feasible as presumably a time at which all available reactive minerals (i.e., which react over appropriate timescales for leakage) are depleted would be reached. The time for depletion would depend on exact sediment composition, which for Vrøgum was not fully confirmed. An increasing magnitude of acidification did however continue until the end of the injection phase concurring with this hypothesis. 5.3. Risks to water resources The reactive transport model showed development of distinct geochemical zones (as described in Section 5.2) each with implications for risks posed to water resources. The initial pulse phase showed the greatest concentrations of dissolved ions and thus infers a risk based on maxima’s achieved. For the Vrøgum aquifer no particularly hazardous trace elements were released meaning risks posed by this phase were small. However, other aquifers may contain minerals with more hazardous constituent elements or surface

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

bound trace elements thus risks will differ between sites. The elevated pulse moved with advective flow inferring it may be the most mobile risk to water resources following leakage of CO2 . Succeeding the pulse phase a zone of aquifer characterized by depletion of the initial buffering mineral (in this case amorphous gibbsite) leads to a zone of increased acidification radiating out from the leakage plane. This zone exhibits reduced dissolved ion concentrations and increasing but still moderate acidification thus would pose little risk to water resources. Finally total exhaustion of buffering capacity is achieved leading to maximal acidification adjacent to the leakage zone. In this zone of advanced acidification other pH sensitive minerals (if present) could now dissolve, or pH dependent desorption occur, potentially mobilizing other harmful trace elements. Field concentrations of Ba and Zn were generally higher than could be explained by the reactive transport model, particularly close to the injection plane. This suggests that the lower pH levels reached in the latter stages of the injection phase (likely caused by buffering exhaustion) may have induced increased dissolution of other unidentified host minerals. For Al concentrations, buffering exhaustion would be in some respects self-limiting. Once gibbsite is depleted and it’s dissolved constituent elements displaced, concentrations would inevitably return toward background levels. Direct evidence of which is unfortunately not provided with the current field results as field injection ceased before this can be clearly assessed.

139

5.4. Depth and calcite As seen in Fig. 10 the effects of CO2 contamination are amplified with depth (i.e., PCO2 > 2 bar). Dissolution of gibbsite is enhanced and subsequently Al driven ion exchange is more severe. This ultimately leads to greater dissolved element maximal in the pulse phase with shorter duration (due to more rapid sediment depletion), in turn giving rise to faster and greater acidification. These results suggest leakage would be easier to detect but also more deleterious to water quality with depth. Addition of calcite bestows significant buffering capacity and high resilience to acidification. For compositions as low as 0.1% calcite, pH is initially limited to 6.2 (Fig. 11a), for which similar effects have been previously reported (Cahill et al., 2013). Buffering capacity is quickly exhausted at such low levels and further acidification achieved under which other pH sensitive reactions may now occur. Larger amounts of calcite (e.g., > 0.5%) provide a greater total buffering capacity, inferring longer resilience to acidification. In general results for calcite inclusion infer leakage may be easier to detect and pose less risk from species mobilized at low pH, typically heavy metals, due to the lower degree of acidification. However this would depend on the solid solution composition of calcite which could contain impurities, e.g., U, which has been observed at leakage natural analog sites as previously described (Keating et al., 2013).

Fig. 12. Schematic representation of likely water chemistry evolution in shallow unconfined silicate dominated aquifer system in space and time during a CO2 leakage scenario (based on observed field and modeling results from Vrøgum). Schematic shows key phases and zones of water chemistry evolution including; pulse phase, post pulse phase, buffering exhaustion interface and advanced acidification with associated generalized chemistry profiles.

140

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141

5.5. MMV program implications Results described here have clear implications for CCS MMV programs. Development of water chemistry is not simple or linear, but rather evolves temporally and spatially (with flow, distance and depth). Current results suggest that for settings like Vrøgum (i.e., silicate dominated) several generalized phases of chemical development should be considered as shown in Fig. 12. At early times the pulse phase would likely be detected first due to advective transport. This phase could be detected most efficiently by increases in alkalinity and elements associated with buffering minerals and exchange processes (e.g., Ca, Mg), i.e., indirectly as electrical conductivity. Within the pulse phase, pH may not be significantly reduced and thus have limited use as a leakage indicator. Following the pulse phase (i.e., in intermediate times) a post-pulse zone forms. This zone may be the most likely detected region of leakage as growth is relatively rapid and sustained, dependent on buffering capacity and groundwater flow. Elevated levels of exchangeable ions are present and may allow confirmation of leakage; however, this would be dependent on aquifer composition and secondary buffering mineral dissolution rate. If the rate of dissolution is too low, aqueous concentrations (and EC) may only be elevated within natural variation thus pH may now become a more reliable indicator of leakage. Following full depletion of buffering capacity, aqueous elemental concentrations (and therefore EC) would return to near background levels. In this zone of the aquifer pH and alkalinity will exhibit the most significant deviation from natural conditions and form the best detection parameters. 6. Conclusions A 1D reactive transport model of physically simulated CO2 leakage was formed in PHREEQC based on comprehensive field measurements. The model was able to simulate water chemistry evolution following injection of gas phase CO2 into a shallow, siliclastic potable aquifer. A range of geochemical processes were included in the model in various combinations with Al driven (gibbsite derived) cation exchange able to simulate the chemistry evolution in time and space most accurately. Alkalinity, pH, Ca, Mg, Na, K, Al, Ba and Sr were all simulated to a relatively accurate degree for the 72 days of continuous measurements provided by the field study. Mn, Fe and Zn were not successfully modeled suggesting other processes are controlling aqueous concentrations, as yet unidentified. The final model was calibrated with finite amounts of amorphous (reactive) and crystalline (less reactive) gibbsite. Glacial sand laboratory batch reactor experiments were also conducted and subsequently modeled. The batch model confirmed the inclusion of gibbsite derived Al driven cation exchange to describe water chemistry change in an independent manner. Following formation of the flow model, future predictions were made in order to demonstrate potential impacts over longer timescales. Model results show water chemistry during CO2 leakage evolves in time and space, thus risks to water resources and leakage detection need to be considered appropriately. In this case study development of zones of acidification were evident whereby depletion of finite amounts of buffering minerals led to increasing acidification radiating out from the leakage plane. This phenomenon poses a risk to water resources due to pH sensitive mineral dissolution and desorption and subsequent increased mobility of pH sensitive toxic trace species. Overall results infer risks from leakage in such a system will change with time, most likely increasing as acidification increases. The most appropriate parameters by which to detect a leak will also change in time with EC more suitable in the early stages and pH in the latter. Increasing depth in the model (inferred as increasing PCO2 ) was seen to amplify observed effects, increasing the speed of sediment depletion and thus acidification. Addition of calcite

bestowed significant buffering capacity, limiting pH decrease (to 6.1) and increasing exchange processes driven by additional Ca in solution. References Assayag, N., Matter, M., Goldberg, D., Agrinier, P., 2009. Water–rock interactions during a CO2 injection field-test: implications on host rock dissolution and alteration effects. Chem. Geol. 265 (1–2), 227–235. Bower, K.M., Zyvoloski, G., 1997. A numerical model for thermo–hydro–mechanical coupling in fractured rock. Int. J. Rock Mech. Min. Sci. 34 (8), 1201–1211. Boyd, A.D., Liu, Y., Stephens, J.C., Wilson, E.J., Pollak, M., Peterson, T.R., Einsiedel, E., Meadowcroft, J., 2013. Controversy in technology innovation: contrasting media and expert risk perceptions of the alleged leakage at the Weyburn carbon dioxide storage demonstration project. Int. J. Greenhouse Gas Control 14, 259–269. Bruggenwert, M.G.M., Kamphorst, A., 1979. Survey of experimental information on cation exchange in soil systems. Dev. Soil Sci. Volume 5 (Part B), 141–203 (Chapter 5). Cahill, A., Jakobsen, R., Mathiesen, T.B., Jensen, C.K., 2013. Risks attributable to water quality changes in shallow potable aquifers from geological carbon sequestration leakage into sediments of variable carbonate content. Int. J. Greenhouse Gas Control 19, 117–125. Cahill, A.G., Jakobsen, R., 2013. Hydro-geochemical impact of CO2 leakage from geological storage on shallow potable aquifers: a field scale pilot experiment. Int. J. Greenhouse Gas Control 19, 678–688. Cahill, A.G., Jakobsen, R., Marker, P.A., 2014. Hydrogeochemical and mineralogical effects of sustained CO2 contamination in a shallow sandy aquifer: a field scale controlled release experiment. Water Resour. Res., 50. Carroll, S., Hao, Y., Aines, R., 2009. Geochemical detection of carbon dioxide in dilute aquifers. Geochem. Trans., 10. Dahmke, A.G., Matthess, A., Pekdeger Schenk, D., Schulz, H.D., 1986. Near-surface geochemical processes in quaternary sediments. J. Geol. Soc. 143, 667–672. Dzombak, D.A., Morel, F.M.M., 1987. Development of a database for modeling adsorption of inorganics on iron and aluminum-oxides. Environ. Prog. 6 (2), 133–137. Engesgaard, P., Jensen, K.H., Molson, J., Frind, E.O., Olsen, H., 11 1996. Large-scale dispersion in a sandy aquifer: simulation of subsurface transport of environmental tritium. Water Resour. Res. 32, 3253–3266. Freyberg, D.L., 1986. A natural gradient experiment on solute transport in a sand aquifer.2. Spatial moments and the advection and dispersion of nonreactive tracers. Water Resour. Res. 22 (13), 2031–2046. Gaus, I., Audigane, P., Andre, L., Lions, J., Jacquemet, N., Dutst, P., Czernichowski-Lauriol, I., Azaroual, M., 2008. Geochemical and solute transport modelling for CO2 storage, what to expect from it? Int. J. Greenhouse Gas Control 2 (4), 605–625. Goldich, S.S., 1938. A study in rock-weathering. J. Geol. 46 (1), 17–58. Hansen, B.K., Postma, D., 1995. Acidification, buffering, and salt effects in the unsaturated zone of a sandy aquifer, Klosterhede, Denmark. Water Resour. Res. 31 (11), 2795–2809. Humez, P., Audigane, P., Lions, J.C., Chiaberge, G., fant Bellenfant, 2011. Modeling of CO2 leakage up through an abandoned well from deep saline aquifer to shallow fresh groundwaters. Transp. Porous Media 90 (1), 153–181. Humez, P.P., Negrel, V., Lagneau, J., Lions, W., Kloppmann, F., Gal, R., Millot, C., Guerrot, C., Flehoc Widory, D., Girard, J.F., 2014. CO2 -water-mineral reactions during CO2 leakage: geochemical and isotopic monitoring of a CO2 injection field test. Chem. Geol. 368, 11–30. Jaffe, P.R., Wang, S., 2003. Potential effect of CO2 releases from deep reservoirs on the quality of fresh-water aquifers. Greenhouse Gas Control Technol., 1657–1660 (Vols I and Ii Proceedings). Jensen, K.H., Bitsch, K., Bjerg, P.L., 1993. Large-scale dispersion experiments in a sandy aquifer in Denmark – observed tracer movements and numerical-analyses. Water Resour. Res. 29 (3), 673–696. Jessen, S.D., Postma, F., Larsen, P.Q., Nhan, L.Q., Hoa, T.K., Pham, T., Long, T.V., Viet, P.H., Jakobsen, R., 2012. Surface complexation modeling of groundwater arsenic mobility: Results of a forced gradient experiment in a Red River flood plain aquifer, Vietnam. Geochim. Cosmochim. Acta 98, 186–201. Keating, E.H., Fessenden, J., Kanjorski, N., Koning, D.J., Pawar, R., 2010. The impact of CO2 on shallow groundwater chemistry: observations at a natural analog site and implications for carbon sequestration. Environ. Earth Sci. 60 (3), 521–536. Keating, E.H., Hakala, J.A., Viswanathan, H., Carey, J.W., Pawar, R., Guthrie, G.D., Fessenden-Rahn, J., 2013. CO2 leakage impacts on shallow groundwater: field-scale reactive-transport simulations informed by observations at a natural analog site. Appl. Geochem. 30, 136–147. Kharaka, Y.K., Thordsen, J.J., Kakouros, E., Ambats, G., Herkelrath, W.N., Beers, S.R., Birkholzer, J.T., Apps, J.A., Spycher, N.F., Zheng, L.E., Trautz, R.C., Rauch, H.W., Gullickson, K.S., 2010. Changes in the chemistry of shallow groundwater related to the injection of CO2 at the ZERT field site, Bozeman, Montana. Environ. Earth Sci. 60 (2), 273–284. Kjøller, C., Larsen, F., Postma, D., 2000. Nickel mobilization in a sandy aquifer in response to groundwater acidification. Groundwater, 259–260. Kjoller, C., Postma, D., Larsen, F., 2004. Groundwater acidification and the mobilization of trace metals in a sandy aquifer. Environ. Sci. Technol. 38 (10), 2829–2835.

A.G. Cahill, R. Jakobsen / International Journal of Greenhouse Gas Control 37 (2015) 127–141 Klusman, R.W., 2011. Comparison of surface and near-surface geochemical methods for detection of gas microseepage from carbon dioxide sequestration. Int. J. Greenhouse Gas Control 5 (6), 1369–1392. Larsen, O., Postma, D., Jakobsen, R., 2006. The reactivity of iron oxides towards reductive dissolution with ascorbic acid in a shallow sandy aquifer – (Romo, Denmark). Geochim. Cosmochim. Acta 70 (19), 4827–4835. Laudelou, H., Vanblade, R., Bolt, G.H., Page, A.L., 1968. Thermodynamics of heterovalent cation exchange reactions in a montmorillonite clay. Trans. Faraday Soc. 64 (546P), 1477. Lemieux, J.M., 2011. Review: the potential impact of underground geological storage of carbon dioxide in deep saline aquifers on shallow groundwater resources. Hydrogeol. J. 19 (4), 757–778. Lions, J., Devau, N., de Lary, L., Dupraz, S., Parmentier, M., Gombert, P., Dictor, M.C., 2014. Potential impacts of leakage from CO2 geological storage on geochemical processes controlling fresh groundwater quality: a review. Int. J. Greenhouse Gas Control 22, 165–175. Lu, J.M., Partin, J.W., Hovorka, S.D., Wong, C., 2010. Potential risks to freshwater resources as a result of leakage from CO(2) geological storage: a batch-reaction experiment. Environ. Earth Sci. 60 (2), 335–348. Navarre-Sitchler, A.K., Maxwell, R.M., Siirila, E.R., Hammond, G.E., Lichtner, P.C., 2013. Elucidating geochemical response of shallow heterogeneous aquifers to CO2 leakage using high-performance computing: implications for monitoring of CO2 sequestration. Adv. Water Resour. 53, 45–55. Parkhurst, D.L., C.A.J. Appelo (2013). Description of input and examples for PHREEQC version 3—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations U.S. Geological Survey Techniques and Methods, book 6, chap. A43. Peter, A., Lamert, H., Beyer, M., Hornbruch, G., Heinrich, B., Schulz, A., Geistlinger, H., Schreiber, B., Dietrich, P., Werban, U., Vogt, C., Richnow, H.H., Grossmann, J., Dahmke, A., 2012. Investigation of the geochemical impact of CO2 on shallow groundwater: design and implementation of a CO2 injection test in Northeast Germany. Environ. Earth Sci. 67 (2), 335–349. Postma, D., Appelo, C.A.J., 2000. Reduction of Mn-oxides by ferrous iron in a flow system: Column experiment and reactive transport modeling. Geochim. Cosmochim. Acta 64 (7), 1237–1247. Postma, D., Larsen, F., Hue, N.T.M., Duc, M.T., Viet, P.H., Nhan, P.Q., Jessen, S., 2007. Arsenic in groundwater of the Red River floodplain, Vietnam: controlling geochemical processes and reactive transport modeling. Geochim. Cosmochim. Acta 71 (21), 5054–5071.

141

Reuss, J.O., Cosby, B.J., Wright, R.F., 1987. Chemical processes governing soil and water acidification. Nature 329 (6134), 27–32. Siirila, E.R., Navarre-Sitchler, A.K., Maxwell, R.M., McCray, J.E., 2012. A quantitative methodology to assess the risks to human health from CO2 leakage into groundwater. Adv. Water Resour. 36, 146–164. Sverdrup, H., Warfvinge, P., 1988. Weathering of primary silicate minerals in the natural soil environment in relation to a chemical-weathering model. Water Air Soil Pollut. 38 (3–4), 387–408. Trautz, R.C., Pugh, J.D., Varadharajan, C., Zheng, L.G., Bianchi, M., Nico, P.S., Spycher, N.F., Newell, D.L., Esposito, R.A., Wu, Y.X., Dafflon, B., Hubbard, S.S., Birkholzer, J.T., 2013. Effect of Dissolved CO2 on a shallow groundwater system: a controlled release field experiment. Environ. Sci. Technol. 47 (1), 298–305. Varadharajan, C., Tinnacher, R.M., Pugh, J.D., Trautz, R.C., Zheng, L.G., Spycher, N.F., Birkholzer, J.T., Castillo-Michel, H., Esposito, R.A., Nico, P.S., 2013. A laboratory study of the initial effects of dissolved carbon dioxide (CO2 ) on metal release from shallow sediments. Int. J. Greenhouse Gas Control 19, 183–211. Viswanathan, H., Dai, Z.X., Lopano, C., Keating, E., Hakala, J.A., Scheckel, K.G., Zheng, L.G., Guthrie, G.D., Pawar, R., 2012. Developing a robust geochemical and reactive transport model to evaluate possible sources of arsenic at the CO2 sequestration natural analog site in Chimayo, New Mexico. Int. J. Greenhouse Gas Control 10, 199–214. Wilkin, R.T., Digiulio, D.C., 2010. Geochemical impacts to groundwater from geologic carbon sequestration: controls on pH and inorganic carbon concentrations from reaction path and kinetic modeling. Environ. Sci. Technol. 44 (12), 4821–4827. Xu, T.F., Spycher, N., Sonnenthal, E., Zhang, G.X., Zheng, L.E., Pruess, K., 2011. TOUGHREACT version 2.0: a simulator for subsurface reactive transport under non-isothermal multiphase flow conditions. Comput. Geosci. 37 (6), 763–774. Zheng, L.G., Apps, N., Spycher, J.T., Birkholzer, Y.K., Thordsen, J., Beers, S.R.W.N., Herkelrath, Kakouros, E., Trautz, R.C., 2012. Geochemical modeling of changes in shallow groundwater chemistry observed during the MSU-ZERT CO2 injection experiment. Int. J. Greenhouse Gas Control 7, 202–217. Zheng, L.G., Apps, J.A., Zhang, Y.Q., Xu, T.F., Birkholzer, J.T., 2009. Reactive transport simulations to study groundwater quality changes in response to CO(2) leakage from deep geological storage. Greenhouse Gas Control Technol. 9 (1(1)), 1887–1894. Zheng, L.G., Spycher, N., Birkholzer, J., Xu, T.F., Apps, J., Kharaka, Y., 2013. On modeling the potential impacts of CO2 sequestration on shallow groundwater: transport of organics and co-injected H2S by supercritical CO2 to shallow aquifers. Int. J. Greenhouse Gas Control 14, 113–127.