On the potential for CO2 mineral storage in continental flood basalts – PHREEQC batch- and 1D diffusion–reaction simulations

Continental flood basalts (CFB) are considered as potential CO2 storage sites because of their high reactivity and abundant divalent metal ions that can potentially trap carbon for geological timescales. Moreover, laterally extensive CFB are found in many place in the world within reasonable distances from major CO2 point emission sources. Based on the mineral and glass composition of the Columbia River Basalt (CRB) we estimated the potential of CFB to store CO2 in secondary carbonates. We simulated the system using kinetic dependent dissolution of primary basalt-minerals (pyroxene, feldspar and glass) and the local equilibrium assumption for secondary phases (weathering products). The simulations were divided into closed-system batch simulations at a constant CO2 pressure of 100 bar with sensitivity studies of temperature and reactive surface area, an evaluation of the reactivity of H2O in scCO2, and finally 1D reactive diffusion simulations giving reactivity at CO2 pressures varying from 0 to 100 bar. Although the uncertainty in reactive surface area and corresponding reaction rates are large, we have estimated the potential for CO2 mineral storage and identified factors that control the maximum extent of carbonation. The simulations showed that formation of carbonates from basalt at 40 C may be limited to the formation of siderite and possibly FeMg carbonates. Calcium was largely consumed by zeolite and oxide instead of forming carbonates. At higher temperatures (60 – 100 C), magnesite is suggested to form together with siderite and ankerite. The maximum potential of CO2 stored as solid carbonates, if CO2 is supplied to the reactions unlimited, is shown to depend on the availability of pore space as the hydration and carbonation reactions increase the solid volume and clog the pore space. For systems such as in the scCO2 phase with limited amount of water, the total carbonation potential is limited by the amount of water present for hydration of basalt.


Introduction
Underground sequestration of carbon dioxide is a potentially viable greenhouse gas mitigation option as it reduces the release rate of CO 2 to the atmosphere [1]. CO 2 can be trapped subsurface by four storage mechanisms: (1) structural and stratigraphic trapping; (2) residual CO 2 trapping; (3) solubility trapping; and (4) mineral trapping [2]. Mineral trapping has been considered as the safest mechanism in long-term storage of CO 2 [3].
Mineral storage of CO 2 in basaltic rocks is favored over siliciclastic reservoirs both by the higher abundance of divalent metal ions in basalt and the faster reactivity of basaltic glass or crystalline basalt [4]. Moreover, basalts such as the Columbia River flood basalts (CRB) are abundant and in many places close to CO 2 point source emissions [5]. During the last decade several flood basalts around the world have been mapped for the possibility of CO 2 storage, and possible candidates such as CRB in USA and the Deccan traps in India have been identified [4][5][6].
To be a candidate for CO 2 storage, the flood basalt must have a proper sealing and sufficient injectivity, the latter limited by the available connected pore space. In flood basalts, the connected pore space is typically found at zones containing abundant vesicles or in breccias between basalt flows. Because central zones of flood basalts commonly are dense and impermeable without vesicles, and flows are laterally continuous over large areas and commonly stacked vertically for hundreds of meters, flow units can act as seals [5]. The non-porous inner parts of flows may however be penetrated by networks of vertical fractures. These fractures can be open and conductive, or closed by mineralization and nonconductive.
The main objectives of this study were to performe batch-and 1D diffusion-reaction numerical simulations to determine the geochemical potential for secondary carbonate formation and to estimate the volume changes and the possibility of self-sealing following the basalt-CO 2 interactions. The CRB system was used as an example case and our results were compared to earlier reported laboratory experiments and numerical simulations of CO 2 -basalt interactions. As CO 2 stored underground will distribute spatially in the reservoir to give a range of reactive conditions, such as the potential of reactions by H 2 O dissolved in supercritical CO 2 [7,8] or reactions in the H 2 O-rich phase from residually trapped CO 2 , we divided the simulations into three systems representing different parts of CO 2 storage: (1) basalt alteration in the H 2 O-rich phase at constant CO 2 pressure; (2) basalt alteration in a H 2 O saturated CO 2 phase, and (3) reactions at the boundary of the CO 2 plume where CO 2 diffuses into the aquifer from the boundary of the CO2 plume ( Figure 1).

Methods
All thermodynamic and kinetic calculations were performed using the geochemical code PHREEQC-2 [9]. This code is capable of simulating complex interactions between dissolved gases, aqueous solutions, and mineral assemblages in batch and 1D advection-diffusion-reaction mode. As the code can only model fully saturated systems, natural systems must be simplified to end-member situations, such as given by constant pressure boundary conditions as may be the case close to underground CO 2 plumes, or the assumption of packages (batches) of water reacting along a reaction path with a homogenous sediment or rock body. Based on these limitations we divided the simulations into three systems representing different parts of CO 2 storage: (1) basalt alteration in the H 2 O-rich phase at constant CO 2 pressure; (2) basalt alteration in a H 2 O saturated scCO 2 phase, and (3) reactions at the boundary of the CO 2 plume where CO 2 diffuses into the aquifer from the boundary of the scCO 2 plume (Figure 1). In the second case, we assumed that the CO 2 phase had swept through the systems and dried out residual water, giving only dissolved water in the scCO 2 phase. In this case an upper limit of carbonation potential was estimated as reactions were allowed to occur until (nearly) all water was consumed, passing the upper 2 mol/Kgw theoretical limit for the Truesdell-Jones activity model [9].
The standard state adopted in this study for the thermodynamic calculations was that of unit activity for pure minerals and H 2 O at any temperature and pressure. For aqueous species other than H 2 O, the standard state was unit activity of the species in a hypothetical 1 molal solution referenced to infinite dilution at any temperature and pressure. For gases, the standard state was for unit fugacity of a hypothetical ideal gas at 1 bar of pressure. All simulations used the llnl.dat database based on the thermo.com. V8.R6.230 dataset prepared at the Lawrence Livermore National Laboratory, with additions of thermodynamic data for those phases not present (see description below). CO 2 fugacity coefficients were estimated according to the modified Redlich-Kwong (SRK) equation of state [10] and the solubility was adjusted for by a poynting correction term (exp(v CO2 (P sat -P)/RT) where v denotes molar volume, P pressure, R the universal gas constant and T absolute temperature) [11]. The density of CO 2 at 40 C and 100 bar was approximated from Bachu and Stewart [12] to be 600 Kg/m 3 and the solubility of water in scCO 2 at the same conditions was approximated to 0.5 mole% [13,14] The simulations were divided into batch simulations of the H 2 O rich and CO 2 rich phases respectively, and 1D diffusion of CO 2 in the H 2 O rich phase to obtain information on the CO 2 -basalt interactions over a continuous range of CO 2 pressures. The latter was solved by PHREEQC using @ t C ¼ D L @ 2 x C þ q , where C denotes molal (mol/Kgw) concentration, q denotes the sink term, subscripts t and x refer to derivatives in time and x-direction respectively, and an efficient diffusion coefficient D L of 0.45x10 -9 m 2 /s was used for CO 2 [15] and all solutes. Sketch of possible reaction settings during CO 2 storage in basalt. System 1 (S1) is close to the injector and contains a wet CO 2 (0.5 mole% H 2 O at 100 bar and 40 C) with no residual water; system 2 (S2) is fully in the H 2 O rich phase with CO 2 diffusing in from the plume boundary; and system 3 (S3) is at the boundary of the CO 2 plume with both sufficient non-wetting CO 2 at a constant CO 2 partial pressure of 100 bar and with sufficient water wetting the mineral surfaces and available for reactions.
Dissolution rates of minerals in the basaltic rock were calculated according to a kinetic equation taking into account pH and the distance from equilibrium: where S is the reactive surface area (m 2 ), k i are rate constants (moles/m 2 s), a H is the H + activity, n is the reaction order with respect to H + and OH -, and Ω is the saturation state given by: Where ΔG r is the Gibbs free energy of the reaction, R is the gas constant, and T is absolute temperature. Reaction rate constants for crystalline basalt (pyroxenes and plagioclase) were obtained from Palandri and Kharaka [16], and pH dependencies were taken from the same source. The dissolution rate of basaltic glass was calculated according to the expression suggested by Gislason and Oelkers [17]: where k + is the far-from-equilibrium dissolution rate coefficient. The saturation state term 1-Ω was approximated to 1 (i.e., rate independent to distance from equilibrium) supported by earlier numerical estimates of glass-CO 2 reactivity suggesting an approximately linear relation between time and reaction progress for basaltic glass [18]. This expression takes into account the effect of pH as well as the effect of the concentration of solutes such as fluoride as they complex with Al 3+ and reduce the Al 3+ activity [19]. The specific surface area for basalt (m 2 /g) was estimated by: where the ratio A p /V p denotes the ratio between pore surface and pore volume (m -1 ), ϕ is connected porosity, and ρ s (g/m 3 ) is the density of the basalt solid estimated from the fraction of the individual basalt components. A S sp value of 1.52 × 10 -5 m 2 /g basalt (= 0.137 m 2 /Kg water) was obtained for the CRB using an average basalt solid density of 2.93 × 10 6 g/m 3 with 10% connected pore space and a A p /V p ratio of 400 m -1 [5]. The reactive surface area was calculated from the mass of the glass and minerals present according to: where M and n are molar mass and moles of mineral i, and X r is the fraction of the total mineral surface that is reactive. As X r is highly uncertain and is suggested to vary by orders of magnitude [20,21], we used a value of 0.1 for the base case and varied X r from 1 to 10 -3 . The use of mass or mass fractions of the individual basalt components to estimate the release rates of elements from the basalt is supported by a recent experimental study which suggests that release rates estimated from the sum of volume fractions of the constituent minerals are within one order of magnitude from measured values [22]. A list of kinetic parameters is given in Table 1. All secondary phases were allowed to form according to the local equilibrium assumption [23]. Changes in solid-phase volumes and porosities ϕ caused by the mineral reactions were calculated according to: where ϕ t=0 is the initial porosity, n and v are moles and molar volume of mineral i respectively, and V total is the total volume of the system. The basalt was defined to consist of a mixture of glass and crystalline basalt with mineral and glass fractions chosen based on reported data from CRB [6,24,25]. To represent the crystalline basalt, plagioclase (Ca 0.5 Na 0.5 Al 1.5 Si 2.5 O 8 ) and the pyroxenes augite (Ca 0.7 Fe 0.6 Mg 0.7 Si 2 O 6 ) and pigeonite (Ca 1.14 Fe 0.64 Mg 0.22 Si 2 O 6 ) were chosen. The hydrolysis equilibrium constants of these phases were estimated using the PHREEQC program assuming ideal solid solutions of the end-members enstatite, ferrosilite and wollastonite for the pyroxenes, and albite and anorthite for the Table 1 Kinetic parameters for dissolution of primary minerals based on empirical data given in Palandri and Kharaka [16] and for basaltic glass from [17]  plagioclase. Equilibrium constants for the solid solutions for temperatures up to 100 C were estimated with PHREEQC and from these data coefficients a to e for the PHREEQC built-in analytical expression (log K = a + bT + c/T + dlog 10 (T) + e/T 2 ) were estimated using nonlinear regression in MATLAB. The glass composition (Ca 0.015 Fe 0.095 Mg 0.065 Na 0.025 K 0.01 Al 0.105 S 0.003 Si 0.5 O 1.35 ) was taken from [6] and modified by adding a small fraction of sulfur which is a common minor constituent of the CR basaltic glass [26].
The secondary mineral assemblage was chosen based on reports on basalt weathering [27][28][29][30], with additional carbonates that could potentially form at elevated CO 2 pressures from the release of Fe, Mg and Ca. The ankerite composition chosen for this work was CaFe 0.6 Mg 0.4 (CO 3 ) 2 which corresponds to a solid solution of 0.6 ankerite (CaFe(CO 3 ) 2 ) and 0.4 dolomite (CaMg(CO 3 ) 2 ). Because ankerite (CaFe 0.6 Mg 0.4 (CO 3 ) 2 ) was not listed in the thermodynamic database, we estimated values using the same approach as in [31]. The full list of secondary minerals is given in Table 2.
To simulate the CRB-CO 2 interaction we used the average concentrations of solutes reported for the Grand Ronde Formation (Table 3). As supercritical CO 2 (scCO 2 at T > 31.1 C; P > 73.9 bar) is the preferred choice for CO 2 storage, based on higher density compared to gaseous CO 2 , we simulated aqueous-phase basalt-CO 2 interaction at a depth of 800 meters at a CO 2 pressure of 100 bar and temperatures of 40 to 100 C. The reactivity of basalt and a H 2 O saturated scCO 2 phase was simulated using an estimated 0.5 mol% H 2 O and a CO 2 density of 600 g/cc giving an initial mass of 0.003 Kg H 2 O per 1 liter pore space.

Results
System 1: Basalt alteration in the H 2 O-rich phase at constant CO 2 pressure i) CRB mineral and glass dissolution and formation of secondary minerals Following the injection of CO 2 into the system, pH immediately decreased from 9.5 to below 4, and thereafter gradually increased to 5.8 at the end of 10000 years (Figure 2a). At the acidic pH secondary phases such as saponite (Ca 0.165 Mg 3 Al 0.33 Si 3.67 O 10 (OH) 2 ), celadonite (KMgAlSi 4 O 10 (OH) 2 ) and zeolite (stilbite) were thermodynamically stable and formed (Figure 2b). Glass dissolved orders of magnitude faster than the crystalline basaltic constituents and more than half dissolved after 10000 years (Figure 2c). The dissolution rate of glass was not increased by the aqueous fluoride as the Al 3+ activity was fixed by the kaolinite and amorphous silica equilibria. The fluoride therefore only increased the total soluble aluminium. The steady release of Fe from the basalt saturated the water with respect to siderite and a total amount of 10 moles/kgw formed after 10000 years ( Figure 2d). Other carbonates, such as ankerite, dolomite, magnesite, and dawsonite did not form as elements such as Mg and Ca was consumed by the non-carbonate secondary phases.
The effect of temperature on the basalt hydration and carbonation was investigated by simulating the system at 60, 80 and 100 C ( Figure 3). As for the 40 C simulation we see that basaltic glass dissolves orders of magnitude faster than the crystalline basalt components and the glass is the major source for the secondary phases. The dissolution rates of the basalt components scale exponentially with temperature, and the glass is almost The mineralogy of the CRB has been described in [25,32] and the weight fraction of pyroxene, feldspar and glass was estimated as average values from the reported data. completely dissolved after 10000 years at 60 C, whereas the time for a complete dissolution takes 4000 and 1500 years at 80 and 100 C respectively (Figure 3a, d, g). The secondary mineral assemblages were largely the same for all temperatures. Stilbite dominated together with amorphous silica (40 and 60 C) and quartz (80 and 100 C) (Figure 3b, e, h). Saponite formed at 40 and 60 C, but not at higher temperatures. Other secondary minerals such as albite, celadonite, and kaolinite formed at all conditions. At 60 C, magnesite and dolomite were still considered to be too slow to form (see [29]) and siderite was the only phase that formed. At 80 and 100 C, magnesite and later ankerite formed together with siderite. Taking zero porosity as the maximum extent of possible reactions we see that the total amount of CO 2 trapped as solid carbonates did not change much with temperature ( Figure 4). The reaction rates increased however with temperature and the time needed to reach the maximum potential therefore decreased with temperature ( Figure 4).
ii) On the limitation of pore-space for the basalt carbonation Secondary phases such as stilbite and amorphous silica have lower density than the basalt components and alteration therefore leads to a reduction of pore space. At the presence of CO 2 , secondary carbonates further reduce the pore space. For the volume calculations we used expression (6) with the molar volumes listed in Table 2. At 40 C, the starting porosity of 10% is reduced to 0.85% after 10000 years. At the higher temperatures all porosity is lost after 2700, 1200, and 300 years respectively at 60, 80, and 100 C ( Figure 5). Taking the extreme of 0% porosity as the limit for the reactions we obtain a maximum carbonation potential (mol CO 2 stored/Kgw) at the different temperatures of 13.5, 29.3, and 28.5 moles for 60, 80, and 100 C ( Figure 5). The simulated clogging of the pore space fits well with shortterm laboratory percolation experiments on open-system   basalt-CO 2 alteration which shows loss of porosity and a rapid reduction of permeability during CO 2 -basalt interactions (e.g., [33]).
iii) Reduction of pore-space as a function of reactive surface area As the reactive surface area is a large uncertainty we simulated the changes of porosity over a range of values from a maximum being equal to the estimated physical surface area S 0 (equation (5) with X r = 1) to a three orders of magnitude reduction ( Figure 6). The physical conditions of the simulated system was the same as for the base-case at 40 C and a CO 2 pressure of 100 bars. At a reactive surface area that is equal to the estimated S 0 all porosity is lost after approximately 1000 years as stilbite and siderite fills the pore space. If the reactive surface area is reduced by one order of magnitude (i.e. the base case) nearly 1/10 of the original 10% porosity is preserved. Further reductions by one and two orders of magnitude lead to smaller changes and at three orders of magnitude reduction relative to S 0 almost no change is observed ( Figure 6).

System 2: The potential for carbonate growth in a H 2 Osaurated scCO 2 phase
The reaction between H 2 O dissolved in scCO 2 and basalt was simulated at 100 bar pressure and 40 C. The initial amount of water was 0.003 Kg and no H 2 O was allowed to enter the system. This is an ideal endmember case and serves to illustrate the carbonation potential in a volume with limited hydration potential.
As secondary phases such as stilbite formed, water was rapidly consumed and most gone after 45 years (Figure 7a). At this point stilbite was unstable and supplied water until all water was consumed after approximately 100 years (Figure 7a). Following the basalt hydration, siderite and ankerite formed from the released Ca, Mg, and Fe, with a final total amount of 0.2 moles CO 2 consumed per liter pore space after 100 years (Figure 7b). If H 2 O had been allowed to dissolve into the scCO 2 phase from residual aqueous phases trapped in the smaller pores, the carbonation potential would have been larger. This process is however, to the knowledge of the authors, not possible to simulate using the PHREEQC code, and was therefore outside the scope of this study.

System 3: 1D diffusion of CO 2 into the CRB aquifer
To see how the basalt reacted under different CO 2 pressures, we defined a 1D diffusion-reaction simulation. This provided us with basalt-CO 2 interactions over a continuous range of CO 2 pressures from the background 1 bar up to the maximum 100 bars. The system corresponds to a stagnant zone presented as a column with one end close to the boundary of the injected CO 2 plume and the other end further away from the plume (Figure 1). The distance reached for the CO 2 into the column is given by the balance between diffusive transport and removal of carbon by secondary carbonate formation. We therefore varied reaction rates from no reaction giving the maximum lengtht of diffusive transport, and up to the base-case rate given by a reaction surface area 1 order of magnitude lower than the estimated physical surface area. Figure 8 shows pH, dissolved CO 2 (mol/Kgw) and amount of secondary carbonate formed in the 1D column. As CO 2 diffuses into the column pH drops to approximately 4 at full  Figure 5). The simulations suggest that the total amount of secondary carbonates that form is dictated by the available pore space and the thermodynamic stability of secondary phases rather than temperature, whereas carbonate generation rates depend on the exponential increase of basalt dissolution rates with temperature. Figure 5 Porosity changes caused by the basalt alteration at 40 to 100 C. Secondary hydrated species and carbonates incorporate the H 2 O and CO 2 masses into the solids and clogs the pore space. As reaction rates increase exponentially with temperature, the pore space is filled up faster at the higher temperatures.  (Figure 8a). The penetration depth decreased rapidly if reactions were allowed as siderite formed and pulled carbon out of aqueous system (Figure 8b, c). At the highest reaction rate (base-case), CO 2 diffused less than 10 meters into the column as approximately 13 moles/Kgw of siderite formed at the end of the 10000 years simulation (Figure 8d). Siderite formed at greater depth if the reactive surface area was reduced by another order of magnitude, but less formed in total (Figure 8d).

Uncertainty on the reactive surface area
The reactive surface area is considered as a major source of uncertainty (e.g., [20,34]) and this leads to corresponding high uncertainties in timing and extent of reactions as dissolution rates have a first order dependence on reactive surface areas. Weathering rates in nature are commonly observed to be 1-3 orders of magnitude lower than in laboratory experiments (e.g., [20,21,34]), and this may in part be explained by differences in reactive and physical (total) surface area between experimental and natural systems. We assumed in this study a base-case reactive surface area 1 order of magnitude lower than the estimated physical surface area for the basalt. A further two orders of magnitude reduction in the reactive surface area, which is within the range of values expected for natural systems, resulted in little basalt alteration and only minor reduction of porosity (see Figure 5). A better understanding of the surface area of porous basalt and the effect of time (aging) on features such as dislocation densities and reactive surface areas are therefore required to understand the potential for CO 2 mineral storage in basaltic rocks.

Uncertainty on the choice of secondary phases used in the model
Growth rate experiments of carbonates such as magnesite and dolomite have shown that the activation energy is high and that growth is negligible at low temperatures (e.g., [35][36][37]). Dissolution rate studies of siderite suggests that the reaction rate is intermediate between calcite and magnesite [38,39], and growth rate data suggest that siderite may form down to room temperature [40]. Data on ankerite dissolution and growth is to the knowledge of the authors not known. The crystallographic and physical characteristics of ankerite do resemble those of dolomite and siderite, and the chemistry is related to dolomite with the Mg 2+ substituted by various amounts of Fe 2+ and Mn 2+ . If the growth rate is close to the magnesian carbonates such as dolomite and magnesite [41,42], the amount that may form during lowtemperature alteration is likely low. In this case, more iron would be available for siderite growth. If on the other hand the growth rate is closer to siderite, we would expect ankerite or other FeMg solid solution carbonates to grow during low-temperature alteration.
One uncertainty related to the local-equilibrium assumption is on the growth retention time for the secondary carbonates. The local-equilibrium assumption predicts growth of the secondary phases as soon as an infinitesimally small supersaturation is reached [23]. The time it takes to nucleate sufficient mass to initiate a significant growth may however be hundreds to thousands of years for some secondary phases [31]. There are no nucleation rate data for siderite and ankerite and the retention time is hence unknown.
Finally, the total potential for secondary carbonate growth may be affected by the amount of magnesium and iron that enters ferromagnesian calcites. As a significant fraction of the metal cations may substitute for calcium (e.g., [43]), a iron-magnesium rich calcite may Figure 7 Basalt alteration in the scCO 2 rich phase with initial 3 grams of water per liter pore space. A) as zeolites form H 2 O is consumed and the water activity is reduced. After approximately 45 years most water is consumed, whereas all is gone before 100 years. B) siderite formed the secondary carbonate initially followed by ankerite.
potentially form rather than ankerite and thereby reduce the amount of siderite formed.

Comparisons to experiments, numerical simulations and natural analogues of basalt-CO 2 interactions
Our simulations suggest that the potential for carbonate growth is limited to siderite or FeMg carbonates at low temperatures as secondary phases such as zeolites outcompeted the carbonates for calcium. We here compare our simulated results with reported data on CO 2 -basalt interactions from laboratory experiments, natural analogues, and other reported numerical simulations.
The reactivity of CRB and other continental flood basalts are available from the long-term (months to years) laboratory experiments done by Schaef and coworkers [6,24]. In these experiments basalt samples from USA, India, South Africa, and Canada were reacted with CO 2 at about 100 bars and 60 to 100 C. Reacted samples from these experiments showed generation of Ca-rich carbonates interpreted as calcites with minor siderite and magnesite. In experiments on CRB using mixtures of H 2 S and CO 2 at 60 C and 100 bar and run for 181 days, pyrite (FeS 2 ) formed together with Mg-Fe poor calcite and a Ca-poor Fe-carbonate [6]. Our simulations at the same temperatures show rapid formation of siderite (60 C) or siderite and magnesite at higher temperatures ( Figure 3). Our simulations do not predict any calcite growth as the calcium activity is lowered by zeolite formation. Calcite would however form in our models if the zeolites were not allowed to form at local equilibrium, and possibly if a magnesian ferroan (solid solution) calcite was used in the model instead of the pure end-member calcite. Therefore, the apparent difference between our model and the experiment may be caused by our use of the local equilibrium assumption, whereas the zeolites in the laboratory experiments did not form at low temperatures due to slow kinetics. Recent experiments on basalt dissolution support the preferential release of Mg and Fe over Ca at acidic conditions [22], suggesting that the MgFe-carbonates will dominate as secondary carbonates during CO 2 storage in basalt.
Our numerical simulations share some similarities to other works such as by Marini [18] and Gysi [44], but Partial pressure of the inlet boundary was fixed at 100 bar and with a column temperature of 40 C. As consumption of CO 2 by siderite growth affects the depth of CO 2 diffusion, we ran a sensitivity study on reactive surface area going from no reaction (a) up to the base case (c). Finally the amount of CO 2 trapped as solid carbonate (siderite) was compared for the base case and reduced reactive surface area (d). We see that the reaction rates strongly constraint the depth of the column affected by the CO 2 diffusion. our model and hence the outcome is different in several aspects. The most comprehensive work done earlier is the numerical simulations done by Marini [18] on the reactivity of crystalline and glassy CFB following CO 2 storage. The initial mineralogy was similar to our study whereas the temperature of 60 C was slightly higher than our base case 40 C. In [18] the CO 2 -basalt interactions were stretched to last for more than 280000 years compared to our 10000 years perspective. The main differences between our model and [18] are on the choice of secondary mineral assemblage, and on the focus of limiting factors such as the availability of water for hydration in the present work. The lack of zeolites and hydrous phases other than kaolinite and goethite in [18] made Ca available for secondary carbonates and the total potential for carbonate formation was higher than in our work. Marini allowed dolomite and magnesite to form at 60 C, whereas our simulations only produced siderite at the similar conditions. Moreover, the formation of dawsonite in [18] is still uncertain and possibly limited at high silica activities and with an assemblage of stable NaAl-silicates defined to form [45]. Based on two different approaches, the reactive surface area for basalt was estimated to quite similar values. We estimated a specific surface area of approximately 1.5x10 -5 m 2 /g basalt (= 0.14 m 2 /Kg water at 10% porosity) based on the A p /V p values estimated by [46] and reported in [5], and reduced this value by one order of magnitude to get the reactive surface area. Marini used a geometric model giving a reactive surface area of 0.41 m 2 /Kg water. The higher reactive surface area and higher temperature of [18] resulted in faster reactions and more rapid clogging of the pore space (within a few years). Studies of natural basalt systems at similar or higher temperatures may give some insight into how fast pore space is clogged by basalt hydration or carbonation, and this should be used to improve the estimates of reactive surface areas of basalt for future studies.
Another numerical study on low-temperature (25 C, 30 bar CO 2 ) basaltic glass alteration was presented by Gysi et al. [44]. Again a main difference is on the choice of secondary minerals. Gysi et al. [44] allowed dolomite, magnesite, and Fe-Mg carbonate to form together with calcite and siderite, whereas we did not allow other Mg-Fe carbonates to form than ankerite. As previously stated, the low-temperature formation of dolomite and magnesite is not likely because of the high apparent activation energy and small kinetic coefficients for the growth of Mg-carbonates [35][36][37]. Other carbonates such as siderite and potentially FeMg-calcites are more likely to form at these low temperatures. The high reactive surface area used in [44] is based on a geometric model for glass fragments, and is hence not directly comparable with the surface area estimated for a vesicle pore space of a solid basalt. Although no inverse modeling was done to estimate the reactive surface area of the basalt in [44], fragmented basaltic rocks such as hyaloclastite breccias are expected to have significantly higher reactive surface areas than porous solid basalts, and they are therefore correspondingly more reactive.
One example of a natural analogue that shed light on CO 2 basalt interactions is the CO 2 charged basalt hosted groundwaters at Hekla, Iceland. Solution aqueous species sampled from natural cold springs and rivers here showed a drop in total inorganic carbon (TIC) that was interpreted to result from considerable formation of secondary carbonate phases such as calcite [47]. Reaction path modeling of the system suggests however that the carbonate formation is associated with high pH in accordance with the low TIC in the sampled waters. This system is therefore different from basalt CO 2 storage projects where higher CO 2 pressures may be maintained over time and the pH is lower. In addition to calcite, dolomite was also suggested as a potential storage host for the low temperature reactions in Hekla [47]. This may however be questionable as long-term laboratory experiments at room temperature have failed to form dolomite even at significant super saturations [48], explained by the high activation energy for dolomite growth [32,41]. Another natural analogue that more closely corresponds to industrial CO 2 storage is the basalt-hosted petroleum reservoir on Nuussuaq, West Greenland. In this system the bulk carbonate formation appears to have occurred as secondary weathering products. Other alteration products such as zeolites and oxides were replaced by dolomite, magnesite, siderite, and calcite at temperatures of 70-120 C [49]. Therefore, taking into account the basalt weathering products and not only primary basalt minerals appears to be vital in estimating the total potential for secondary carbonate formation and the long-term potential for CO 2 storage in basalt systems.

Summary and conclusions
Simulations of closed-system (P CO2 = 100 bar, 40 C) and 1D reaction-diffusion (P CO2 = 0-100 bar, 40 C) alteration of basalt suggest that the potential of secondary carbonate formation is limited to siderite at low temperatures as divalent metal cations are preferentially consumed by zeolites and oxides. Higher temperatures 60 -100 C appear to be in favor of secondary carbonate formation, allowing the precipitation of carbonates such as magnesite, siderite and possibly dolomite and other FeMg carbonates (ankerite). Given an unlimited source of CO 2 (fixed CO 2 pressure), the total amount of CO 2 stored as solid carbonates is orders of magnitude higher than the 1-2 mol/Kg water solubility in the formation water ( Figure 4). The total amount trapped might however be reduced if CO 2 , H 2 O or pore space are limiting factors. The formation of secondary hydrous and carbonate phases increases the volume of solids and the porosity is correspondingly reduced ( Figure 5). This together with the immobilization of CO 2 by solid carbonate formation is in favor of safe long-term storage of CO 2 in basaltic aquifers.

Competing interests
We have no competing interests with any organization to publish the manuscript 'On the potential for CO2 mineral storage in continental flood basalts -PHREEQC batch-and 1D diffusion-reaction simulations'.
Authors' contributions VTH Pham has made substantial contributions to conception and designs the manuscript. Her contributions are to acquisition of data, analysis and interpretation of data. PA did revising and has given final approval of the version to be published. HH involved in drafting the manuscript, writing methodology part and revising of the final version. All authors read and approved the final manuscript.