Multidimensional reactive transport modeling of CO2 mineral sequestration in basalts at the Hellisheidi geothermal field, Iceland

Multidimensional reactive transport modeling of CO 2 mineral sequestration in basalts at the Hellisheidi geothermal field, Iceland E.S.P. Aradottir a,b,∗ , E.L. Sonnenthal c , G. Bjornsson d , H. Jonsson b a Reykjavik Energy, Baejarhalsi 1, 110 Reykjavik, Iceland b Science Institute, University of Iceland, Dunhaga 3, 107 Reykjavik, Iceland c Lawrence Berkeley National Laboratory, 1 Cyclotron Rd, Berkeley, CA 94720, USA. (Sonnenthal is partially supported by the U.S. Department of Energy and LBNL under Contract No. DE-AC02-05CH11231.) d Reykjavik Geothermal, Sudurlandsbraut 18, 108 Reykjavik, Iceland Two and three-dimensional field scale reservoir models of CO 2 mineral sequestration in basalts were developed and calibrated against a large set of field data. Resulting principal hydrological properties are lateral and vertical intrinsic permeabilities of 300 and 1700 × 10 −15 m 2 , respectively, effective matrix porosity of 8.5% and a 25 m/year estimate for regional groundwater flow velocity. Reactive chemistry was coupled to calibrated models and predictive mass transport and reactive trans- port simulations carried out for both a 1200-tonnes pilot CO 2 injection and a full-scale 400,000-tonnes CO 2 injection scenario. Reactive transport simulations of the pilot injection predict 100% CO 2 mineral capture within 10 years and cumulative fixation per unit surface area of 5000 tonnes/km. 2 Correspond- ing values for the full-scale scenario are 80% CO 2 mineral capture after 100 years and cumulative fixation of 35,000 tonnes/km 2 . CO 2 sequestration rate is predicted to range between 1200 and 22,000 tonnes/year in both scenarios. The predictive value of mass transport simulations was found to be considerably lower than that of reactive transport simulations. Results from three-dimensional simulations were also in significantly better agreement with field observations than equivalent two-dimensional results. Despite only being indicative, it is concluded from this study that fresh basalts may comprise ideal geological CO 2 storage formations. 1. Introduction Reducing anthropogenic CO 2 in the atmosphere is among the great challenges of this century (e.g. Broecker, 2007; Hoffert et al., 2002; Lackner, 2003; Pacala and Socolow, 2004; Oelkers and Schott, 2005, and references therein). Large scale capture and seques- tration of atmospheric CO 2 has been proposed in an attempt to solve this problem and in situ mineral sequestration is among the suggested methods (Metz et al., 2005; Schrag, 2007). In in situ mineral carbonation, injected CO 2 is combined with dissolved met- als to form carbonate minerals, such as calcite (CaCO 3 ), dolomite (CaMg(CO 3 ) 2 ), magnesite (MgCO 3 ) and siderite (FeCO 3 ), that are stable and environmentally benign over geological time scales. For- mation of these minerals can enhance the stability of injected CO 2 substantially (Metz et al., 2005; Oelkers et al., 2008). Great advances have been made in geochemical computer simulations in recent years and decades. Numerous geochemical ∗ Corresponding author at: Reykjavik Energy, Baejarhalsi 1, 110 Reykjavik, Iceland. Tel.: +354 516 6992; fax: +354 516 6709. E-mail address: edda.sif.aradottir@or.is (E.S.P. Aradottir). Computer codes have been developed, enabling scientists to peer underground and study the complex interplay of dissolved constituents and fluids by other means than are common in the laboratory. Initially, geochemical simulations involved time- independent equilibrium speciation of multi-component systems (e.g. Wolery, 1979; Reed, 1982; Arnorsson et al., 1982). Reaction path models that capture the sequence of chemical/mineralogical states followed the time-independent equilibrium calculations, providing a way for quantitatively and sequentially interpreting irreversible processes such as chemical weathering and hydrother- mal alteration (e.g. Helgeson, 1968; Helgeson et al., 1969). Lichtner (1988) linked mass transport and travel time to the reaction path approach. This approach has now been generalized to het- erogeneous multi-dimensional and multi-component flow fields. Resulting reactive transport models can be applied to fundamental Earth Sciences research, as they provide an integrated way to look at the big picture in the complex interplay of material flow, trans- port and reactions that characterizes subsurface processes, along with assessing the relative importance of individual sub-processes (Steefel et al., 2005). Such complex interplay is to be expected in association with CO 2 mineral sequestration in basaltic rocks and reactive transport models should therefore be ideal for simulating


Introduction
Reducing anthropogenic CO 2 in the atmosphere is among the great challenges of this century (e.g. Broecker, 2007;Hoffert et al., 2002;Lackner, 2003;Pacala and Socolow, 2004;Oelkers and Schott, 2005, and references therein). Large scale capture and sequestration of atmospheric CO 2 has been proposed in an attempt to solve this problem and in situ mineral sequestration is among the suggested methods (Metz et al., 2005;Schrag, 2007). In in situ mineral carbonation, injected CO 2 is combined with dissolved metals to form carbonate minerals, such as calcite (CaCO 3 ), dolomite (CaMg(CO 3 ) 2 ), magnesite (MgCO 3 ) and siderite (FeCO 3 ), that are stable and environmentally benign over geological time scales. Formation of these minerals can enhance the stability of injected CO 2 substantially (Metz et al., 2005;Oelkers et al., 2008).
Great advances have been made in geochemical computer simulations in recent years and decades. Numerous geochemical Computer codes have been developed, enabling scientists to peer underground and study the complex interplay of dissolved constituents and fluids by other means than are common in the laboratory. Initially, geochemical simulations involved timeindependent equilibrium speciation of multi-component systems (e.g. Wolery, 1979;Reed, 1982;Arnórsson et al., 1982). Reaction path models that capture the sequence of chemical/mineralogical states followed the time-independent equilibrium calculations, providing a way for quantitatively and sequentially interpreting irreversible processes such as chemical weathering and hydrothermal alteration (e.g. Helgeson, 1968;Helgeson et al., 1969). Lichtner (1988) linked mass transport and travel time to the reaction path approach. This approach has now been generalized to heterogeneous multi-dimensional and multi-component flow fields.
Resulting reactive transport models can be applied to fundamental Earth Sciences research, as they provide an integrated way to look at the big picture in the complex interplay of material flow, transport and reactions that characterizes subsurface processes, along with assessing the relative importance of individual sub-processes . Such complex interplay is to be expected in association with CO 2 mineral sequestration in basaltic rocks and reactive transport models should therefore be ideal for simulating Fresh 25 • C groundwater, used for dissolving CO2, will be pumped from well HN-01 and transported to well HN-02 in a pipeline (blue line). CO2 captured at the pilot gas separation plant will be transported in a pipeline (yellow line) to well HN-02 where it will be injected along with water from HN-01. Deviation of wells HN-04 and HK-34 is depicted by red lines. Combed lines represent faults. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) its effect in the subsurface. In addition, the models provide tools to predict and optimize long-term management of injection sites, to quantify the amount of CO 2 that can be mineralized, and identify secondary minerals that compete with carbonates for cations leached from the rock matrix.
In this article, we present development of two and threedimensional field scale reactive transport models that simulate mass and reactive transport, along with the mineral alteration associated with injecting dissolved CO 2 into basalts. The models were built using the TOUGH2 (Pruess, 1991), iTOUGH2 (Finsterle, 1999) and TOUGHREACT (Xu et al., 2006(Xu et al., , 2011 simulators. These simulators have been extensively used in simulations of processes associated with geothermal systems and geological storage of CO 2 .

The CarbFix project
CarbFix (e.g. Aradóttir et al., 2011;Gíslason et al., 2010;Matter et al., 2008;Oelkers et al., 2008) is a combined industrial/academic collaboration project between Reykjavík Energy, the Institute of Earth Science at the University of Iceland, Earth Institute-Lamont-Doherty Earth Observatory at Columbia University in New York and the Centre National de la Recherche Scientifique/Universite Paul Sabatier in Toulouse, that was developed in order to assess the feasibility of in situ CO 2 mineral sequestration in basaltic rocks in Iceland. The project consists of a CO 2 pilot injection, laboratory based experiments, study of natural analogs, predictive model development, numerical modeling and model validation, as well as cost analysis. Rather than injecting CO 2 directly into geological formations, CarbFix is developing a technology to dissolve CO 2 into formation fluids and well water during injection. Once dissolved, CO 2 is no longer buoyant compared to pore fluids, improving considerably security due to decreased leakage risks. This approach of solubility trapping also promotes carbonation of the host rock and thus facilitates the safe long-term sequestration of CO 2 in the subsurface. Fig. 1 shows an aerial map of the CarbFix injection site. The site is located ca. 3 km SW of the Hellisheidi geothermal power plant, owned and operated by Reykjavík Energy. The power plant is currently generating 303 MWe and 133 MWth by tapping water and steam from a high temperature reservoir at 800-3000 m depth. The power plant annually captures some 40,000 tonnes CO 2 of volcanic origin as a by-product of the geothermal energy production. Reykjavík Energy has already invested substantially in a comprehensive network of deep and shallow wells in the area, in piping and pumping systems and in a pilot plant that separates CO 2 from other volcanic gases.
The wellfield shown in Fig. 1 comprises ideal conditions for studying the long-term effect of mineral CO 2 sequestration in basaltic rocks. Connection to Hellisheidi power plant provides access to a concentrated source of CO 2 . Additionally, field data, such as drill cuttings and a calcite rich cap-rock overlying the high-temperature reservoir, suggests that mineral CO 2 sequestration already plays an important role in the evolution of the Hellisheidi geothermal system (see, e.g. Níelsson and Franzson, 2010;Gebrehiwot et al., 2010). Injecting and precipitating CO 2 into nearby formations with the objective of imitating and accelerating this natural CO 2 sequestration process should therefore be considered as an environmentally benign process. Wiese et al. (2008) carried out a study of CO 2 fixation in high-temperature geothermal systems in Iceland and found natural CO 2 sequestration in Hellisheidi to be in the range of 4100-23,500 tonnes/year or about 650,000 tonnes/km 2 . The study showed sequestration to be most intense from 550 to 800 m below surface.

Process description
A pilot gas separation station has been built next to Hellisheidi geothermal power plant. The pilot station will separate CO 2 from the geothermal gas coming from the condensers of the power plant by sequential extraction (Aradóttir et al., 2011). In the pilot study, approximately 5% of the total geothermal gas coming from the power station will be separated, resulting in the capture of approximately 2000 tonnes CO 2 for CarbFix on an annual basis. After its capture, CO 2 will be transported in a 3 km pipeline to the pilot injection site as high pressurized gas. There, CO 2 will be injected at depth into well HN-02. Co-injected water will divert injected CO 2 further downward, resulting in a single phase fluid entering the 400-800 m targeted storage formation (Aradóttir et al., 2011).
CO 2 injection was started in Hellisheidi in January 2012. Environmental authorities in Iceland have granted licenses for the injection. It is anticipated that the results of the CarbFix project will be used to optimize the in situ carbon mineralization process, enabling it to also work at other basalt rich sites worldwide.

Injection site infrastructure
The injection site in Hellisheidi is outfitted with ten wells, whose location is shown in Figs. 1 and 3. Fresh groundwater, used for dissolving CO 2 , will be pumped from well HN-01 and transported to well HN-02 in a 1 km pipeline. There it will be injected along with CO 2 from the pilot gas separation plant. The main permeable section in HN-02 is located close to 500 m depth and most of the injected fluid should enter the reservoir there. The static water level in HN-02 is at about 100 m depth and its injectivity by gravity only is about 50 kg/s. Wells HK-12, HK-25, HK-07 and HK-13 serve as monitoring wells and only penetrate the 5-10 • C, shallow, free surface and fresh groundwater system in the topmost 200 m. Wells HN-04, HK-34, HK-31 and HK-26, however, are drilled down through the storage formation and serve as deep monitoring wells. This reservoir is acting as confined in the short term and extensive well-logging indicates it to follow an approximately linear temperature gradient of 100 • C/km (see, e.g. Björnsson et al., 2006;Gunnarsson et al., 2011;Franzson et al., 2004). Wells HN-04, HK-34, HK-31 and HK-26 are all situated downstream from the injection well. Wells HK-31 and HK-26 are located about 1750 and 3000 m downstream from well HN-02, while the main permeable zones in the directionally drilled HN-04 and HK-34 are at about 70 and 300 m lateral distance from injection well HN-02, respectively. Well HN-04 is most permeable around 400 m depth, while the main permeable section of well HK-34 is located slightly below 500 m depth.

Modeling studies within CarbFix
Numerical modeling plays an important role in the CarbFix project as it provides tools to predict and optimize long-term management of the injection site as well as to quantify the amount of CO 2 that can be mineralized. A series of models have been developed within CarbFix, as summarized in Table 1. TOUGH2 (Pruess, 1991), iTOUGH2 (Finsterle, 1999) and TOUGHREACT (Xu et al., 2006(Xu et al., , 2011 have mainly been used for developing the field scale models presented in this article. The workflow of the model development was divided into the following five interrelated stages: (1) development of a conceptual model of regional hydrological flow conditions and reactive transport with respect to ongoing water-rock interaction, (2) compilation and validation of thermodynamic and kinetic databases, (3) parameters in hydrological model constrained by forward and inverse calibration to field data, (4) forward predictive model simulations with a fully coupled reactive mass transport model, and (5) model validation by comparison to laboratory and/or field data. The workflow stages are correlated to the individual model studies in Table 1. Aradóttir et al. (2011) and Aradóttir et al. (in press) discuss in detail the conceptual model of the pilot CO 2 injection, and the compilation and validation of a thermodynamic database suitable for this study. This article, therefore, only gives a brief overview of stages (1) and (2) in the model development, while focusing on stages (3)-(5).

Overview of conceptual and simple reservoir models
Figs. 2 and 3 show, respectively, in a cross-section, a schematic representation of the CarbFix CO 2 injection site and the conceptual reservoir model of the sequestration site (model II in Table 1). Fresh flood basalts govern the bedrock from surface down to 200 m depth. These basalts host a cold (5-10 • C) groundwater system, with the static water table located about 100 m below surface. Below lies a 200 m thick hyaloclastite layer that separates the cold groundwater system from the storage formation below due to significantly lower vertical permeability. The storage formation extends from about 400 to 800 m depth. It is composed of relatively fresh basaltic lavas interbedded with minor hyaloclastite layers/intrusions. Highest permeability is encountered around 500 m depth.
Two slug type tracer tests have been carried out at the injection site since 2007 with the objective of studying regional hydrology in the area and seeing if tracer dispersivity could identify reservoir flow type (fracture or matrix) within the storage formation. Fig. 4 shows tracer breakthrough curves in well HN-04 from both tracer tests. Rezvani Khalilabad et al. (2008) used a one dimensional and multipath tracer dispersion model to study recovery of the first tracer test (model I in Table 1). The model predicted most of the storage formation to consist of a large volume of relatively homogeneous porous media. A low volume and fast breakthrough path was found to channel only about 3% of the tracer flow between wells HN-02 and HN-04, and thus only plays a minor role in the storage formation.
When the primary rock comes in contact with the injected weakly acidic, carbonated fluid, dissolution reactions occur, leaching cations out of the rock matrix. As dissolution reactions proceed, protons are consumed, the pH of formation fluids increases, and concentration of leached cations builds up. At certain concentrations, the water becomes supersaturated with respect to secondary minerals like carbonates, which begin to precipitate. The efficiency of this reactive chemical transport is highly sensitive to the surface area available for water-rock interaction.
Natural analogs provide important insight into the secondary mineralogy associated with CO 2 rich and CO 2 depleted water-basalt interaction and hence, which minerals are likely to compete with carbonates for dissolved cations in in situ CO 2 mineral sequestration (model III in Table 1). Aradóttir et al. (2011) cover in detail the alteration mineralogy associated with basalt weathering under low and elevated CO 2 conditions. Simple oxides and hydroxides, clay minerals, carbonates and zeolites are among common alteration minerals.

Flow conditions
In this study, two and three-dimensional models of the Carb-Fix injection site were developed (models V, VI, VII, VIII and IX in Table 1). The models are aligned in the direction of regional groundwater flow, which is from north-northeast to south-southwest. Figs. 3 and 5 show layer distribution in the three-dimensional model and the geographical location of the CarbFix numerical models, respectively. The three-dimensional model consists of eight layers, from 100 to 900 m depth, whereas the two-dimensional model is a simplified version of the three-dimensional model, consisting of a single horizontal layer between 425 and 525 m depth. Each layer consists of 6198 elements. Apart from varying thickness, all layers have identical element geometry. Generally, layers containing feed zones in the three-dimensional model are set narrow but those that do not contain feed zones are thicker. Elements at the north-northeast and south-southwest boundaries are set at constant temperature, pressure and chemical composition, as well as the top and base layers in the three-dimensional model. The mesh is rectangular, except around injection well HN-02, where it is radial in order to get reasonable distribution of injected fluids and dissolved constituents in the near vicinity of the well. Elements close to the injection well are small (<1 m) but their size gradually increases further away from the well. The horizontal hydraulic gradient that drives natural groundwater flow through the bedrock has been estimated to be 5 m/1000 m, making pressure at the models north-northeast boundary 3.36 bar higher than at the south-southwest boundary. Estimation of the hydraulic gradient is based on water table measurements in wells HN-01, HN-02, HN-04, HK-31 and HK-26 that extend over several years. Stratigraphy of the CarbFix model is based on the conceptual model shown in Fig. 3. In the topmost 800 m, rock formations mainly consist of relatively fresh basaltic hyaloclastite formations and lava flows. Hydrological parameters of the lavas were calibrated using field data, but porosity and permeability values for the hyaloclastite formations were taken from Frolova et al. (2005), and set to be 10% and 30 mD, respectively.
Reaction induced matrix permeability changes are calculated from changes in porosity using ratios of permeabilities calculated from the Carman-Kozeny relation (Bear, 1988), ignoring changes in grain size, tortuosity and specific surface area. Aradóttir et al. (2011) carried out an extensive literature review of water-basalt interaction at low and elevated CO 2 conditions with the objective of predicting which minerals are likely to dissolve and precipitate in in situ CO 2 mineral sequestration. Subsequently, Aradóttir et al. (in press), developed and validated a thermodynamic database describing mineral reactions of interest for CO 2 -water-basalt interaction. All primary and secondary minerals in the database were added to the reactive transport models and their dissolution reactions written in terms of the same basis species as used by Aradóttir et al. (in press) and Gysi and Stefánsson (2011), i.e. H 2 O, H + , HCO − 3 , Al(OH) − 4 , Fe 2+ , Mg 2+ , Na + , K + , SiO 2(aq) ,

Fig. 2.
A simplified schematic representation of the CarbFix pilot CO2 injection in Hellisheidi, Iceland. Regional groundwater flow will carry injected carbonated water downstream. Cations leached from the rock matrix by the weakly acidic carbonated plume will react with CO 2(aq) and form carbonates.
Cl − and O 2(aq) . The minerals are listed in Table 2. Selection of basis species is primarily based on the findings of Stefánsson et al. (2001), who studied dominant aqueous species in natural waters in the basaltic terrain of Iceland. Initial water composition of all model elements at time zero is that of well HN-01 in Hellisheidi, which is given in Table 3. Carbonated injection water has the same composition, except that CO 2 concentration correspond to about Henry's law saturation at 30 bar CO 2 partial pressure, resulting in approximately 1 M total concentration of H 2 CO 3 , HCO − 3 and CO 2− 3 . As a result of acidification, pH of the boundary water drops down to 3.6. Initial rock is composed of plagioclase, pyroxene, olivine and basaltic glass. Basaltic glass composition is that of Oelkers and Gíslason (2001), whereas the specific composition of plagioclase, pyroxene and olivine are those given by Stefánsson (2001). Table 4 shows the initial mineral volume composition used in model simulations.

Kinetics of mineral dissolution and precipitation
Dissolution and precipitation of all minerals except allophane, Al(OH) 3(am) , p-antigorite, calcite, Fe(II) and Fe(III) hydroxides, imogolite and kaolinite is kinetically controlled. Kinetic rates are a product of the rate constant and reactive surface area, according to the following general rate expression, which is transition state theory based Steefel and Lasaga, 1994): where r is rate of dissolution or precipitation, k is the temperature dependent rate constant, A is the specific reactive surface area, K is the equilibrium constant for the dissolution/precipitation reaction taking place and Q is the reaction quotient. Â and Á must be determined by experiment but are often set to unity. As Fig. 3. The CarbFix conceptual model, drawn along A-A in Fig. 1. CO2 is dissolved in water during injection into HN-02, resulting in a single phase fluid entering the storage formation (400-800 m depth). Black arrows show direction of regional groundwater flow and red stars main feed-zones in wells. Layers in the three dimensional CarbFix model are designated with horizontal black lines. Layers that contain feed zones are narrower than others. Pressure at the north-northeast boundary is 3.36 bar higher due to the horizontal hydraulic gradient that drives natural groundwater flow through the bedrock. Temperature and pressure gradients are based on extensive well-logging, which show a near-linear increase of temperature and pressure from surface and down to the top (cap-rock) of the high temperature geothermal system.

Table 2
Chemical composition of the minerals and aqueous species considered in this study. See Aradóttir et al. (2011) and Aradóttir et al. (in press) for further information on mineral selection and thermodynamic data. dissolution and precipitation of minerals are often catalyzed by H + (acid mechanism) or OH − (base mechanism), the rate constant in Eq.
(1) is often written as the sum of three mechanisms: where nu, H and OH denote neutral, acid and base mechanisms, respectively. E A is activation energy, k 25 the rate constant at 25 • C, R is the gas constant, T absolute temperature and a activity of a species. Gíslason and Oelkers (2003) measured and summarized previously measured basaltic glass dissolution kinetics far from equilibrium at temperature ranging from 6 to 300 • C and pH between 1 and 11. The resulting rate law describing basalt glass dissolution is: The rate law was implemented in TOUGHREACT using a general form of a species dependent rate constant that is coded in TOUGHREACT: Initial rock mineral composition and potential alteration minerals along with parameters for calculating kinetic rates for mineral dissolution and precipitation.  Gíslason and Oelkers (2003). Rate constant was renormalized with respect to BET surface area and reduced by 3 orders of magnitude. d From Palandri and Kharaka (2004). e From Gíslason et al. (1997). f From Rimstidt and Barnes (1980). Precipitation occurs under the free energy rate law of Carroll et al. (1998). g From Knauss et al. (2005). h All zeolites assumed to have the same rate law as heulandite from Palandri and Kharaka (2004).
where i denotes the species dependent mechanism and j specific species to which the rate constant depends on. Parameters used for the kinetic rate expression of different minerals are given in Table 4. Forsterite, albite and diopside dissolution data from Palandri and Kharaka (2004), based mostly on data from Pokrovsky and Schott (2000), Chou and Wollast (1985) and Knauss et al. (1993), were used for olivine, plagioclase and pyroxene, respectively. Rate constants for albite were multiplied by a factor of five to account for the effect of anorthite in the plagioclase solid solution as anorthite and albite have very different kinetic behaviors. The same does not apply to olivine and pyroxene end-members. The rate law of the most abundant end-member was therefore chosen to be representative of olivine and pyroxene solid solutions. The rate law of Gíslason and Oelkers (2003) was used for basaltic glass but the geometric rate constant reported by Gislason and Oelkers (see Eq. (3)) was renormalized with respect to BET surface area.
Rate constants of all primary minerals and glasses were reduced by 3 orders of magnitude, as simulations with the reported rate constants gave way too rapid dissolution during preconfiguration of model VI in Table 1. Simulations of rain water interacting with the primary rock composition given in Table 4 were carried out with different rate reduction factors and resulting water composition compared to the chemical composition of groundwater of a known age after corresponding simulation times. Simulations using 3 orders of magnitude reduction factor gave excellent match with field observations and the same rate reduction factor was therefore implemented to the rate constants of the primary minerals in Table 4. An equivalent way of reducing simulation rates would have been to reduce mineral reactive surface areas by three orders of magnitude. Use of a rate reduction factor of this sort is justifiable in light of the fact that experimentally determined weathering rates are commonly two to four orders of magnitude faster than field-derived rates as has, e.g. been shown by White and Brantley (2003), White et al. (1996) and Brantley (1992).
Rate law parameters for moganite and quartz were taken from Gíslason et al. (1997) but from Rimstidt and Barnes (1980) for SiO 2(am) . Rate law parameters for siderite are those reported by Knauss et al. (2005) but rate law parameters for other minerals are from Palandri and Kharaka (2004) who summarized previous rate measurements from the literature and refitted the data. All zeolites were assumed to have the same rate law as heulandite.
Mineral reactive surface areas in the subsurface are generally unknown. In the current study, reactive surface area of primary minerals and glasses was assumed to be 20 cm 2 /g, which is in agreement with the work of Sonnenthal et al. (2005), who calculated reactive surface areas of several minerals assuming a cubic array of truncated spheres constituting the rock framework. Resulting values were roughly in the range of 10-100 cm 2 /g. Surface areas of precipitating minerals are generally much higher than those of the primary rock. In the current study, surface areas of secondary clay minerals, zeolites and carbonates were assumed to be 10,000, 1000 and 500 cm 2 /g, respectively. Secondary SiO 2(s) minerals were assumed to have a surface area of 1000 cm 2 /g.
When the aqueous phase supersaturates with respect to a certain secondary mineral, a small volume fraction of 1 × 10 −6 is used for calculating a seed surface area for the new phase to grow (Xu et al., 2010). The precipitation of secondary minerals is represented using the same kinetic expression as that for dissolution, except for SiO 2(am) which precipitates under the free energy rate law of Carroll et al. (1998). As precipitation rate data for most minerals are unavailable, parameters for neutral pH rates were employed to describe precipitation. This is a critical but necessary assumption because of lack of data on precipitation kinetics.

Simulations
Simulations carried out in this study can roughly be divided into two categories. First, is pure hydrological and tracer recovery modeling aimed at determining primary reservoir properties such as porosity and permeability (model VI in Table 1). Substantial volume of initial state and transient field data was available for constraining the model parameters. Reactive chemistry was then coupled to the calibrated hydrological model and forward predictive simulations carried out on CO 2 -water-basalt interaction associated with pilot and full scale CO 2 injection at Hellisheidi (models VII, VIII and IX).
The numerical models were allowed to reach natural state flow by running no injection/production simulations for 10,000-100,000 years or until pressure and temperature changes between time-steps became negligible. Steady state output was used as initial conditions for other simulations. Batch geochemical simulations of water-rock interaction (model IV) were also carried out before starting fully coupled reactive mass transport simulations in order to pre-react the measured water composition of well HN-01 to a steady state with the primary minerals listed in Table 4.

Model calibration and well production scheme design
Hydrological properties of the storage formation were calibrated using iTOUGH2 (Finsterle, 1999) to simulate tracer tests that have been ongoing in Hellisheidi since 2007 (model V). The iTOUGH2 simulator is well suited for such simulations as it is capable of solving fluid and heat flow both in forward and inverse mode. The calibration was built on seven datasets on water level fluctuations due to sudden injection and/or production stops as well as downstream tracer breakthrough and recovery curves. The EOS1 equation of state was applied along with the 'two water' option of TOUGH2. This approach was chosen as it involves the simplest way of adding a conservative tracer to a flow problem. The resulting model is based on an integrated form of the mass and energy conservation equations (Lichtner, 1996;Björnsson et al., 2006).

Predictive mass transport and reactive transport simulations
Reactive chemistry was coupled to the calibrated hydrological model and TOUGHREACT used for predictive reactive transport simulations. The TOUGHREACT simulator solves the same fluid and heat flow equations as iTOUGH2 but also solves equations describing chemical reactions. Simulations were carried out in two (models VII and IX) and three dimensions (model VIII) and involved the pilot injection test as well as a full-scale injection of all CO 2 currently emitted from Hellisheidi geothermal power plant.
In simulations of the pilot CO 2 injection, 3 kg/s of carbonated water (p CO 2 ∼ 30 bar) were modeled to be injected continuously into well HN-02 for 6 months. Well HN-04 was modeled to produce 1 kg/s continuously for the first 12 months of the pilot test, whereas 1 kg/s production out of well HK-34 was turned on 6 months into the pilot test and went continuously on for 18 months (see Fig. 7). In simulations of a full-scale CO 2 sequestration, 40 kg/s of carbonated water (p CO 2 ∼ 30 bar) were injected continuously into HN-02 for 10 years. No downstream wells were produced in this model. Well HN-01 provided water for injection in both modeling scenarios. CO 2 concentration of injected water was set to be at Henrys law saturation at indicated CO 2 pressure as all gas injected into well HN-02 dissolves during injection.
Mass transport simulations, equivalent to the reactive transport simulations described above, were carried out in order to identify the difference between aqueous concentrations in mass transport and reactive transport simulations. The mass transport simulations only involve mass and heat transport but do not take chemical reactions into account. Dissolved CO 2 is expected to be retained by solid phases in CO 2 mineral sequestration, causing downstream concentrations to be lower in reactive transport simulations. Comparing results of equivalent mass transport and reactive transport simulations provides an assessment of how reliable predictive mass transport simulations are for a reactive system. Such information is important because fully coupled reactive transport simulations on large, complex systems can be impractical to run due to heavy computational requirements. In such cases, mass transport simulations have to suffice.

Model calibration and well production scheme design
Calibration of the hydrological field model (V) with field data resulted in horizontal permeability of lava flows to be 300 mD and a vertical permeability of 1700 mD (mD is milli Darcy, 1 × 10 −15 m 2 ). Active matrix porosity was estimated to be 8.5%. The fact that basalt flows in Iceland are commonly composed of vertical basalt columns surrounded by vertical fracture networks could explain why higher vertical permeability was obtained in calibration of hydrological parameters. One way to incorporate this heterogeneity into the numerical model involves using an idealized dual-porosity model that considers fractures in a porous medium (Xu and Pruess, 2001). The method was, however, not used in the current study due to its heavy computational requirements. Fig. 6 shows comparison between measured and simulated values for two of the calibration datasets used in inverse simulations. Match between measured and simulated values is satisfactory considering that the single porosity model used cannot account for heterogeneities due to, e.g. fractures providing fast flowing but low volume paths between injection and observation wells.
As calibrated hydrological parameters of the storage formation predict groundwater velocity to be low (25 m/year), plans call for accelerating groundwater flow during the pilot CO 2 injection by producing downstream wells. The calibrated numerical model is a valuable tool for designing optimal injection and production schemes (model VI). Simulations showed that production out of HN-04 significantly speeds up downstream recovery of nonreactive tracers, whereas producing HK-34 does not have the same effect. Therefore, additional expenses due to production out of HK-34 are hardly justifiable until injected dissolved constituents have reached well HN-04. Darcy's law is fully applicable due to the low groundwater flow in the storage formation. Its Reynolds number is orders of magnitude lower than 1 (Fetter, 2001).

Predictive mass transport and reactive transport simulations
The current version of the TOUGHREACT code only allows for single central processing unit (CPU) simulations, which causes difficulties in field scale reactive transport simulations due to the computational requirements of such simulations. Three dimensional mass transport and reactive transport simulations of the 10 year long full-scale CO 2 injection turned out to be too computationally intensive for the single CPU version of TOUGHREACT to handle, as the continuous 40 kg/s injection of carbonated water required small time steps for reasonable distribution of injected species. Three-dimensional reactive transport simulations of the pilot CO 2 injection took several months to run, and were thus only barely practical. Fig. 7 summarizes predicted HCO − 3 concentrations (total dissolved DIC) in wells HN-02, HN-04, HK-34, HK-31, HK-26, and HN-01 for both the pilot and the full-scale CO 2 injections. Kinks in curves belonging to the pilot injection are due to start of production out of HK-34 after 6 months and stop of production out of HN-04 after 1 year. In the case of the full-scale injection, injected CO 2 reaches water supply well HN-01 after 4 years of continuous production. Fig. 7 will be discussed in detail in the subsequent sections.

Pilot injection -3D
In the case of the pilot injection, three-dimensional reactive transport simulations predict the CO 2 contaminated plume to reach well HK-34 after 3-4 years. No plume effects are observed in other observation wells. Fig. 7 shows that only a slight difference exists between HCO − 3 concentrations predicted by mass transport and reactive transport simulations in the near vicinity of the injection well (HN-02), but that the difference increases as the CO 2 plume travels downstream. This indicates CO 2 to be retained by solid phases. The predictive value of mass transport simulations of a reactive system therefore decreases as plumes are carried downstream. Dimensionality effects were studied by comparing HCO − 3 concentrations predicted by equivalent two and three-dimensional simulations. Fig. 7 exhibits that a considerable difference exists between concentrations due to dimensionality effects. The threedimensional field model predicts much earlier plume arrival times in downstream wells than the two-dimensional one. Comparison with measured tracer breakthrough curves shows that good agreement exists between measured tracer arrival times and the plume arrival times predicted by the three-dimensional model, which hence are more realistic than those predicted by the twodimensional model. Vertical permeability, which was calibrated to be five times higher than the horizontal permeability, is not accounted for in the two-dimensional model, as it consists of a single horizontal slice. This is likely the main reason for slower breakthrough in the two-dimensional model. Fig. 8 shows cross sectional view of reservoir HCO − 3 concentration after 1 year, 5 years, and 10 years in the CO 2 pilot injection, as predicted by three-dimensional reactive transport simulations. The figure shows that the CO 2 plume travels faster deeper in the reservoir (below 500 m depth) than at shallower depths. This is, e.g. due to production out of HN-01 which supplies solvent water for CO 2 injected into HN-02. Increased diffusion due to higher random kinetic energy at higher reservoir temperature could also play a role. Plume concentration decreases after CO 2 injection is stopped due to dilution, diffusion and chemical reactions, resulting in the plume to be disappearing after 10 years. Fig. 9 shows primary and secondary mineral abundance predicted by three-dimensional reactive transport simulation of the pilot CO 2 injection along with estimated reservoir sequestering efficiency (model VIII). Initially, the primary rock constituted of 43 vol%, 39 vol%, 13 vol%, and 5 vol% plagioclase, pyroxene, olivine and basaltic glass, respectively, of the reservoir (ca. 6.5 km × 4 km × 0.8 km). The volume percentages correspond to 1.0 × 10 13 , 1.5 × 10 13 , 7.2 × 10 12 , and 2.8 × 10 12 moles, respectively. The reservoir bedrock is relatively fresh, but not free of secondary minerals. Clays, simple hydroxides, calcite, quartz, and zeolites are among commonly found alteration minerals in well cuttings collected in the reservoir (e.g. Alfredsson et al., 2007;Helgadóttir et al., 2009). Secondary minerals were not included in the primary rock composition, as quantitative data was not available. Initial interaction between the primary rock minerals and pore water does however, lead to a 'correction' of the primary rock composition toward a more inherent state, which includes clay minerals, simple hydroxides and oxyhydroxides, calcite and SiO 2 phases (SiO 2(am) , quartz and moganite). The fact that the primary rock evolves to its inherent state, suggests both the model and its thermodynamic database to be of proper quality.
Initially, all primary phases dissolve and hydroxides, carbonates and silicon oxides precipitate as the system evolves to its inherent state. Plagioclase and basaltic glass dissolution gradually increases throughout the simulation time, whereas dissolution of pyroxene and olivine levels out. After 5 years, plagioclase dissolution has become considerably greater than dissolution of other primary phases.
The abundance of simple hydroxides and SiO 2 phases (SiO 2(am) , quartz and moganite) increases steadily throughout the whole simulation. Zeolites begin to form after 1 month and rapidly become the second most abundant alteration mineral group in the reservoir. Carbonate formation increases gradually and all the 1100 tonnes of CO 2 injected are predicted to have been sequestered after 10 years. Calcite is most predominant among carbonates but dolomite and a Fe-rich magnesite-siderite solid solution also form to a lesser degree.
Calcite precipitation is most abundant below 500 m depth for the first 5 years but then it spreads out toward 400 m depth. Highest calcite concentrations are observed between wells HN-04 and HK-34. Dolomite precipitates in large quantities at the front of traveling CO 2 plume but dissolves again as the plume passes by. Magnesite-siderite solid solution precipitates within the CO 2 plume. Calcite and magnesite-siderite phases remain stable thereafter. p-Antigorite (Gunnarsson et al., 2005) precipitation is generally associated with calcite precipitation in the simulations. Hydroxides form equally through lava flows in the reservoir and aluminum hydroxide is most abundant. Increased silica precipitation is observed behind the moving CO 2 plume. Amorphous silica is most abundant but quartz and moganite also form. The horizontal zeolite zones so commonly found in Iceland and initially described by Walker (1960) are reproduced in the simulations. Mesolite and heulandite are the most abundant zeolites. Fig. 10 shows a cross sectional reservoir view of selected precipitated minerals in the pilot injection 5 years from injection, as predicted by three-dimensional reactive transport simulations (model VIII). Fig. 11 shows primary and secondary mineral abundance predicted by two-dimensional reactive transport simulation of a 10 year long full-scale CO 2 injection of 40,000 tonnes/year (model IX). Initial interaction between the primary rock and pore water result in the primary rock to evolve to its inherent state, as happened in the pilot simulation discussed above. Results from two-dimensional simulations of the full scale injection somewhat differ the three-dimensional ones from the pilot injection, e.g. due to different CO 2 abundance in the reservoir concentrations. Temperature effects most likely also play a role, as the two-dimensional model only consists of a single slice between 400 and 500 m depth, and thus involves lower temperature than are found at greater depth in the three-dimensional model. This could, e.g. explain the higher fractional abundance of clay minerals in two-dimensional simulations than in three-dimensional runs. Celadonite and other clay minerals are typically more common above 500 m depth in the CarbFix reservoir. Another possible explanation for higher fractional clay abundance in two-dimensional simulations compared to three-dimensional simulations is the initial dominant plagioclase dissolution, which is likely to result in high Al activity and clay formation.

Full scale injection -2D
Basaltic glass dissolution is not observed until after 1 year. Delay of glass dissolution compared to three-dimensional simulations are due to lower simulation temperature in the two-dimensional layer. Plagioclase dissolution gradually becomes most extensive, as was also observed in the three-dimensional pilot injection simulations. Clay minerals are abundant and zeolite precipitation becomes extensive after approximately 10 years. Clays are initially Mg-rich (Mg-chlorite and celadonite), but their Fe-content gradually increases until clay solid solutions become equally Fe-and Mg-rich. Carbonates precipitation involves calcite, dolomite and magnesite-siderite solid solution. Calcite and dolomite are the most and least abundant carbonates, respectively. About 80,000 tonnes CO 2 have been sequestered into carbonates after 100 years. As the two-dimensional model only considers a conservative sequestration volume consisting of a 100 m thick slice between 400 and 500 m depth, it is likely to underestimate the sequestering capacity of the reservoir. The storage formation is four times thicker, so optimally it can be assumed to sequester four times more CO 2 or 320,000 tonnes CO 2 in 100 years.

Discussion
The reactive transport simulations carried out in this study indicate in situ CO 2 mineral sequestration in basalt to be among viable solutions to manage anthropogenic CO 2 emissions. The models predict precipitation of alteration minerals to occur at a considerable distance from the injection well, which indicates that clogging in its near vicinity should not be a problem. CO 2 sequestration rate is modeled to range between 1200 and 22,000 tonnes/year, a value in good agreement with the natural sequestration rate of the Hellisheidi geothermal system, which Wiese et al. (2008) estimated to be 4100-23,500 tonnes/year. Furthermore, models predict cumulative mass of CO 2 fixed in the storage formation per unit surface area in the pilot and full scale injections to be approximately 5000 tonnes/km 2 and 35,000 tonnes/km 2 , respectively. Natural fixation per unit storage area in Hellisheidi was estimated by Wiese et al. (2008) to have accumulated to 650,000 tonnes/km 2 over the past 70,000-140,000 years. The sequestering capacity of the storage formation in Hellisheidi is therefore likely to become significantly greater with time than predicted by 100 year long simulations.
Anthropogenic CO 2 emissions in 2010 were projected to be 35 Gtonnes by Friedlingstein et al. (2010). In global context the results of the current study thus indicate that on an annual basis some 50,000-100,000 km 2 of a few hundred meter thick basaltic formations are needed in order to sequester all anthropogenic CO 2 emissions. When put in perspective to the very large volumes of dry land basalts accessible worldwide, it is not unrealistic to presume that mineral sequestration of CO 2 can develop into being one of the industrially viable solutions available in mitigating the atmospheric CO 2 levels. However, as sequestration was found to be most efficient between approximately 40-80 • C, injection wells should be drilled to a corresponding depth. The storage formation in Hellisheidi also is able to disperse the carbonated water injected, thus maximizing surface area of the reactive transport. Such favorable conditions are harder to find as basalts get buried and older.
Simulated alterations are in agreement with observations from natural analogs. Initial alteration minerals include those commonly found in the CarbFix reservoir at Hellisheidi, which are typical for water-basalt interaction at low CO 2 conditions. Carbonates that form are predominantly calcite -the only carbonate commonly found in Icelandic basalts. In reservoir subdomains containing high CO 2 concentration and low pH, a Fe-rich magnesite-siderite solid solution also forms along with dolomite at the plume front, and these results are in agreement with findings from the basalt hosted petroleum reservoir in Greenland (Rogers et al., 2006). Clay minerals, simple oxides and oxyhydroxides are predicted to compete with magnesite-siderite solid solution for Mg and Fe leached from primary minerals, whereas zeolites compete with calcite for dissolved Ca.
Figs. 9 and 11 show that plagioclase dissolves in largest quantities of primary minerals, and in fact its dissolution becomes orders of magnitude greater than for other minerals with increasing simulation time. Recent work by Gudbrandsson et al. (2011) shows that at low pH, olivine dissolution dominates plagioclase dissolution, which was not the case in the lower temperature two-dimensional simulations carried out in the current study. Therefore, indications are that modeled plagioclase dissolution is too extensive, and as a result thermodynamic and kinetic relations for plagioclase need to be re-evaluated. The excessive plagioclase dissolution leads to a rise in pH, and a build-up of Al, Si, Ca, and Na ions within the aqueous phase. As a result, Al-hydroxide, zeolites and potentially also clays, become supersaturated and precipitate in what is likely to be too massive amounts.
Extensive monitoring during and after the CO 2 pilot injection will provide the basis for the testing, validating and further calibrating the two-and three-dimensional mass transport and reactive transport models discussed in this contribution. Plans call for drilling an exploratory well into the storage formation after the pilot CO 2 injection. The results of reactive transport simulations suggest that it would be wise to wait for about 3-5 years to drill such a well as then collected core samples will be more likely to contain larger amounts of sequestered carbonates and associated alteration minerals. Furthermore the simulations indicate that the exploratory well should be drilled in the north-eastern vicinity of the main feed zone of well HK-34.
The current study indicates how reactive transport models can be used for deepening understanding of subsurface processes associated with CO 2 injection and sequestration as well as designing management schemes for injection sites. However, it is clear that if reactive transport simulations of complex subsurface processes are to become more practical, a parallel version of codes such as TOUGHREACT is essential. Such tools are also necessary if complex calibration of, e.g. reactive surface areas, rate-law parameters, and stochastic permeability is to become feasible.

Uncertainties associated with reactive transport modeling
Large uncertainties are generally associated with reactive transport modeling, in part due to uncertainties in laboratory measured values, but also due to significant differences between values measured in the field and those measured in the laboratory. Often one must thus accept significant uncertainties in parameters describing permeability, diffusivity, reactive surface area and mineral dissolution/precipitation rates. Hence, it follows that the results of reactive transport calculations can have uncertainties as high as several orders of magnitude. Much of this uncertainty can, however, be overcome by obtaining extensive and specific system and/or site specific physical and chemical parameters, as was done in the current study.
A better understanding of reaction mechanisms, and in particular precipitation mechanisms, is an important factor in decreasing uncertainties associated with reactive transport modeling. Discrepancies between rates measured in the laboratory and those measured in the field also need to be further identified. Information on reactive surface areas in the field is also scarce, and often difficult to measure, as it changes with time due to chemical reactions. Xu et al. (2005) performed sensitivity simulations regarding ankerite and dawsonite precipitation kinetics, which indicated changes in reaction rates to result in small changes in precipitated amounts. The reason for this is that precipitation requires reactants whose availability is controlled by the slow dissolution of aluminosilicate minerals.

Conclusions
A series of modeling studies has been carried out with the objective of quantitatively studying CO 2 mineral sequestration in basalts. Two and three-dimensional field scale reservoir models of a proposed CO 2 sequestration formation at Hellisheidi geothermal field were developed through a multistage numerical modeling approach. The models were calibrated against a large set of field and laboratory data. Resulting principal hydrological properties are lateral and vertical intrinsic permeabilities of 300 and 1700 × 10 −15 m 2 , respectively, effective matrix porosity of 8.5% and a 25 m/year estimate for the regional groundwater flow velocity.
Predictive numerical modeling results guide in decision making, indicating, e.g. that it is feasible to pump downstream observation wells in order to expedite recovery of the injected carbonated water. Thereby, the planned pilot CO 2 injection study will generate urgently needed field data to validate reactive numerical models. Predictive reactive transport simulations for both a 1200-tonnes pilot CO 2 injection and a full-scale 400,000-tonnes CO 2 injection scenario show that injected carbonated water is capable of leaching cations from the basaltic rock matrix. This eventually results in pore water becoming super-saturated with respect to carbonates, into which CO 2 is sequestered. Reactive transport simulations of the pilot injection predict 100% mineral capture of injected CO 2 within 10 years, whereas two-dimensional simulations of the full-scale injection predict as much as 80% capture within 100 years.
The predictive value of mass transport simulations was found to be considerably lower than that of equivalent reactive transport simulations as injected CO 2 is retained by solid phases. Furthermore, results from three-dimensional simulations were in significantly better agreement with field observations than results from equivalent two-dimensional simulations. Simulated alterations are in agreement with observations from natural analogs, providing confidence in the developed models and some verification of their accuracy. Clay minerals, simple oxides and oxyhydroxides are predicted to compete with magnesite-siderite solid solution for Mg and Fe leached from primary minerals, whereas zeolites compete with calcite for dissolved Ca.
Despite being only indicative it is concluded from this study that relatively fresh and glassy basaltic formation of high porosity may comprise ideal conditions for geological storage of CO 2 .

DISCLAIMER
This document was prepared as an account of work sponsored by the United States Government. While this document is believed to contain correct information, neither the United States Government nor any agency thereof, nor The Regents of the University of California, nor any of their employees, makes any warranty, express or implied, or assumes any legal responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by its trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof, or The Regents of the University of California. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof or The Regents of the University of California.
Ernest Orlando Lawrence Berkeley National Laboratory is an equal opportunity employer.