Numerical simulations were carried out for determining the chemical reactions relevant for the sequestration of CO2 in basaltic rock formations. The mineralogy of natural geological systems consists of silicate minerals such as the phyllosilicates and zeolites that form complex solid solutions. Using the GEMS code package based on Gibbs energy minimization, combined with the new MINES database, we can now simulate the solubility of these multicomponent and multisite mineral solid solutions in basaltic rocks. This study explores the varying effects of CO2 partial pressures, basaltic glass dissolution kinetics and reaction time on the complex chemistry of the overall CO2-water-basalt reaction path. The simulations indicate four reaction progress stages with the competing reactions between smectites (di- and trioctahedral) and Ca-Fe-Mg-carbonates controlling the amount of CO2 mineralized. A better understanding of these key mineral-fluid reactions and improvement of their thermodynamic models is critical for making more acurate predictive calculations. This comprises the basis for extending the simulations to reactive transport models, and for the assessment of the feasibility of long-term CO2 storage in basaltic rock formations.
Rapid anthropogenic CO2 emissions have raised concerns on climate change, and several carbon capture storage options have been proposed, including in situ trapping in basaltic rock formations , . This geochemical method, also called mineral carbonation, involves the dissolution of the primary basalt forming minerals (i.e. olivine, pyroxene and plagioclase) and basaltic glass followed by the mineralization of CO2 into secondary carbonates. Mineral carbonation is considered as a possible solution for permanently trapping CO2 due to the rapid dissolution rates of basaltic rock forming minerals and the availability of large volumes of basaltic lavas such as the Columbia River basalts in the U.S. and the Deccan Traps in India , . The Carbfix project , in SW Iceland, is an example of a promising field project, which aims at injecting and storing CO2 from a geothermal power plant into the surrounding basaltic rock formations. This field study indicates that 95% of the injected CO2 mineralized within less than 2 years .
Future engineering methods, aiming at optimizing the mineralization of CO2 in basalts, will require a detailed knowledge of the controlling mineral-water reactions and their effects on chemical mass fluxes over different geological timescales. The low temperature alteration mineralogy of a typical CO2-poor natural basaltic system (e.g. Teigarhorn, Eastern Iceland) consists of amorphous and poorly crystalline Fe- and Al-hydroxides, Mg–Fe-bearing smectites and Al–Si minerals (allophane and/or imogolite), and/or chalcedony, Ca–Na zeolites and calcite . An alteration mineralogy consisting of quartz, Fe–Mg carbonate solid solutions and dolomite was found as a result of the interaction of CO2-rich water with basaltic rocks in the Disko-Nuussuaq Island, West Greenland , . Laboratory basaltic glass solubility experiments successfully reproduced this mineralogy at 75°C and pCO2 of ~25 bar within 124 days, whereas at 40°C the mineralization included chemically similar but less crystalline phases , . These experiments showed that the alteration mineralogy and chemistry of the aqueous solutions were largely dependent on acid supply (pCO2), temperature, fluid/rock ratios and reaction time.
We are still facing many challenges to better understand the long-term storage of CO2 in basaltic rock formations. Numerical simulations can be used to evaluate processes occurring at equilibrium vs. the ones driven by kinetics , , , for example, by comparing results to laboratory experiments and natural systems. In a previous numerical modeling study , we used the law of mass action (LMA) code PHREEQC  to simulate chemical reactions associated to CO2-water-basalt interaction reaction path. The simulations indicated that the competing reactions between complex smectite compositions and Fe–Mg-carbonates at low pH (~3–6), and zeolites, smectites and calcite at high pH (~8–10), controlled the availability of Ca, Mg and Fe in aqueous solutions for mineral carbonation. The study of  has focused on simulating the physical aspects of CO2-water-basalt interaction using TOUGHREACT , which included multidimensional reactive transport models to simulate the CO2 sequestration mechanisms in relation to the CarbFix field scale reservoir study mentioned above. One of the main limitations resides in the simulation of more complex non-ideal chemical systems including multicomponent mineral solid solutions where cations can mix on multiple crystallographic sites. Many numerical modeling studies use the LMA method and compile thermodynamic data for minerals in the form of equilibrium constants (logK), while having to use mineral solid solutions with fixed compositions for simulating the alteration of basalts , , . In the present study, I use the GEMS code package (http://gems.web.psi.ch), based on Gibbs energy minimization (GEM), and present the new open access MINES database (http://tdb.mines.edu) for simulating the chemistry of such complex rock alteration reactions.
Some of the challenges that we will review for optimizing the mineralization of CO2 into basalts include adequate mineral solid solutions models, mineral dissolution kinetics and available thermodynamic data for aqueous species and minerals (carbonates, smectites and zeolites). A better understanding of these key mineral-water reactions and improvement of their thermodynamic models is critical for making predictive calculations at the field scale. This comprises the basis for extending future simulations using more complex mineral precipitation/dissolution kinetics, and for the implementation of reactive mass transport models  and pore scale processes  in the assessment of the feasibility of long-term carbonate mineralization in basaltic rock formations.
Reaction path model
Gaseous CO2 dissolves in natural waters to form carbonic acid, which ionizes to bicarbonate at mildly acidic pH conditions according to,
where H2CO3* represents the species CO2(aq) and minor H2CO3 (aq). Reaction of the proton with the primary basaltic minerals (i.e. olivine, pyroxene and plagioclase) and basaltic glass will release divalent cations (i.e. Ca2+, Mg2+ and Fe2+) for the mineralization of CO2 into secondary minerals according to,
Reaction (4) shows that the precipitation of calcite will depend on pH (i.e. increasing pH will drive the reaction to the right). The major problem results from the formation of other secondary minerals than carbonates, which may compete for the divalent cations, especially zeolites that may compete for Ca2+, and the phyllosilicates that may compete for Ca2+, Mg2+ and Fe2+. This series of reactions may be studied in terms of a water-rock reaction path. Reaction path simulations are equivalent to a titration model where the irreversible dissolution of primary rock forming minerals is followed by the precipitation/dissolution of secondary minerals. These reactions can be described by a series of partial equilibria that are similar to the reaction path model of feldspar hydrolysis , . This overall rock dissolution process leads to a pH increase similar to an acid-base titration , where the basalt will act as the base by consuming protons and releasing cations upon dissolution, and the water with its acids such as CO2 will act as the acid. The discrete addition of rock is also called reaction progress in this titration model, where at each step the reaction progress variable is described by addition of an amount of rock, dissolved solutes and a changing secondary mineralogy that is characteristic for the extent of reaction , .
Basaltic glass model and dissolution kinetics
A basaltic rock comprises a rock matrix consisting of volcanic glass and phenocrysts of olivine, pyroxene and plagioclase, and minor Fe-Ti-oxides. Experimental studies have shown that volcanic glasses dissolve much faster than primary rock forming minerals , , and are therefore ideal for prospective CO2 sequestration projects. Since basaltic glass has a bulk composition that is similar to crystalline basalts, this study will use basaltic glass as a proxy for basaltic rocks. This approach is similar to the studies of , , and avoids the complexity of including mineral dissolution kinetics for multiple primary mineral types, which can be added at a later stage once the chemical reactions in the simplified system are better understood. The far-from-equilibrium geometric surface area normalized basaltic glass dissolution rate (r+,geo) was determined experimentally by ,
where AA is a pre-exponential factor including the rate constant (10−5.6 mol of Si/cm2/s), EA an activation energy of 25.5 kJ/mol, R the ideal gas constant, T temperature in K, and ai the activity of the ith aqueous species. The overall forward dissolution rate (r+) based on the geometric reactive surface area (sgeo,t in cm2) at time t (s) can be described by,
The reactive surface area is difficult to estimate in natural geologic rock formations. A simplified assumption is to simulate the reactive surface area assuming spherical particles according to,
where sgeo is the initial reactive surface area in cm2, mt the mole amount at time t and m0 the initial mole amount of basaltic glass, respectively.
Many complexities arise for the formulation of a rate law for the dissolution of this type of multicomponent phase because it will initially dissolve incongruently, i.e. release some cations such as Ca2+ and Na+, whereas higher charged cations such as Si4+ and Al3+ will form stronger bonds within the glass structure . This has led to several theories explaining the dissolution of basaltic glass, including the formation of a hydrated leached layer on the basaltic glass surface through which other cations will diffuse , , , . Approach to equilibrium with respect to the basaltic glass (or its leached layer) and the aqueous solution, will also affect the far-from-equilibrium dissolution rate. This can be evaluated by using a chemical affinity term (1–Q/K), which is a simplified form based on transition state theory , where Q is the reaction quotient and K the equilibrium constant for the solubility of the basaltic glass.
Methods and numerical model
Reaction path modeling
Numerical simulations of the CO2-water-basalt reaction paths were carried out in the thermodynamic system Na–K–Ca–Mg–Fe–Si–C–O–H using the GEMS code package . These simulations are different than in our previous study , where we used the LMA method. In the present study, I make use of the GEM method for solving chemical equilibria, which permits to include temperature and pressure corrections for thermodynamic properties of minerals, gases and aqueous solutions, includes equations of states for gases, and allows calculating non-ideal and ideal multicomponent mineral solid solutions with varying compositions. The starting basalt composition was represented by a natural basaltic glass from Stapafell (SW Iceland), which was used in several previous experimental studies , , , , . The simplified chemical formula of this natural basaltic glass can be described by: K0.008Na0.08Ca0.27Mg0.26Fe(II)0.169Fe(III)0.012SiAl0.35O3.286. The selected starting aqueous solution had a composition of natural spring water (Vellankatla, SW Iceland) buffered by basaltic rocks, which was also used in the experiments of , ,  and the LMA simulations of , . The composition of this aqueous solution is: 256 μmol/kg Si, 269 μmol/kg Na, 11.9 μmol/kg K, 71 μmol/kg Ca, 38 μmol/kg Mg, 0.16 μmol/kg Fe, 1.09 μmol/kg Al, 120 μmol/kg Cl . Its initial composition was varied by saturating it at 50°C with different initial CO2 gas partial pressures (pCO2) ranging between 20 and 100 bar.
Two different simulation types were carried out at 50°C (Table 1): (i) a basalt titration model, and (ii) a kinetic model. In model (i), 0–400 g of basaltic glass was added to the CO2-saturated aqueous solution, where it was assumed that at each step all phases are at equilibrium and that the system undergoes a series of partial equilibria, representing the reaction path of the CO2-water-basalt interaction. In model (ii), 0–150 g of basaltic glass was allowed to dissolve according to the dissolution rate described by Equation (5). The initial geometric reactive surface area was taken as 250 cm2/g , which corresponds to spherical particles with a radius of 45–125 μm. In both of these models, the effects of changing pCO2 and fluid/rock ratio on the reaction path were simulated. In addition, the effects of changing the reactive surface area and time were also simulated in the kinetic model (ii).
|Titration model (i)||Kinetic model (ii)|
|pCO2 (bar)||20, 40, 100||20, 40, 100|
|Initial mCO2 (mol/kg)||0.351, 0.624, 1.056||0.352, 0.628, 1.068|
|Initial pH||3.53, 3.36, 3.20||3.92, 3.68, 3.46|
|Basaltic glass (g)||0–400||150|
|Initial specific sgeo (cm2/g)||–||250|
|EA (activation energy) (kJ/mol)||–||25.5|
|AA (rate constant) (mol Si/cm2/s)||–||10−5.6|
Thermodynamic dataset and activity model
The thermodynamic dataset used in this study is the new MINES database (http://tdb.mine.edu), which consists of a project file that can be managed with the GEMS code package . A detailed list of minerals, aqueous species and gases used in the simulations can be found in the Supplementary Materials Tables 1 and 2. This dataset comprises mineral data of , , and aqueous species from SUPCRT92 (slop98.dat; , , , , ), and was updated with the internally consistent dataset of  for K, Na and Si aqueous species optimized using the GEMSFITS code package , and with more recent experimental data for Al aqueous species , , . The properties of zeolites were derived from the study of , and the properties of phyllosilicates (i.e. saponite, beidellite and nontronite) using the approach from the study of  as described below. The standard state adopted in this study was unit activity for pure H2O and for the aqueous species in a hypothetical one molal standard solution referenced at infinite dilution at any pressure and temperature. The standard state for minerals was that of a pure phase. All activity models and equations of state were calculated using the TSolMod library class implemented in GEM-Selektor . The activity model used for aqueous species was the extended Debye-Hückel according to  and using the Debye-Hückel parameters of , . The equation of state used for real gases (i.e. CO2–CH4–CO–H2O) was the Peng-Robinson-Stryjek-Vera (PRSV) model , .
The thermodynamic properties of di- and trioctahedral smectite end members were estimated using the combined approach of , based on oxide summation techniques and the relation between entropy and volume of compounds. The trioctahedral smectites consisted of the Ca–, Fe–, Mg–, K– and Na– end member varieties of saponite, and dioctahedral smectites of the Ca–, K–, Mg–, Na– end member varieties of nontronite and beidellite (Supplementary Materials Table 1). The standard Gibbs energy was estimated for each of these end members according to the method of , the entropy according to , and the heat capacity function according to . The activities of mineral solid solutions (i.e. carbonates, zeolites and smectites) were approximated using an ideal mixing model, which accounts for changes in the configurational entropy resulting from random mixing of cations on different crystallographic sites (e.g. interlayer and octahedral sites for smectites). One of the advantages of the GEMS code package over our previous simulations study , is that these solid solution models can be defined in the thermodynamic database, and simulations do therefore not necessitate a priori knowledge of the solid solution compositions as is the case for LMA type codes. For the titration model, it was necessary to determine the standard thermodynamic properties of the basaltic glass. The approach taken in this study, was to retrieve them using the basaltic glass formula normalized to one Si and the above mentioned oxide summation techniques and entropy-volume relations used for the phyllosilicates , , .
Results and discussion
Results of the simulations of basaltic glass reacted at 50°C with the CO2-rich aqueous solutions saturated initially at pCO2 of 20, 40 and 100 bar, respectively, are shown in Figs. 1 and 2 . The simulated pH increased upon reaction with the basaltic glass, and displayed multiple stages at which the pH was buffered by several different mineral assemblages. Initially, the pH was buffered at a value of ~3.5 by the ionization of CO2 (aq) leading to the formation of bicarbonate (HCO3−) and proton release (stage 1 in Fig. 1). Upon increased dissolution of the basaltic glass (stage 2 in Fig. 1), the pH reached values of ~5–5.5 and was buffered by a mineral assemblage consisting of chalcedony (SiO2), kaolinite, and Ca–Mg–Fe carbonates (dolomite, Mg–Fe carbonate solid solutions and ankerite). Further basaltic glass dissolution (stage 3 in Fig. 1), led to pH values of ~8–8.5 and a system buffered by the mineral assemblage chalcedony (SiO2), kaolinite, smectites (di- and trioctahedral), and Ca–(Fe)-carbonates (i.e. calcite and ankerite). Stage 3 was reached in all three simulations including the aqueous solutions initially saturated with a pCO2 of 20, 40 and 100 bar. However, a final pH value of ~9.5–10 (stage 4 in Fig. 1) was only reached in the simulation with an initial pCO2 of 20 bar. At these conditions, the stable mineral assemblage consisted of chalcedony (SiO2), kaolinite, smectites (di- and trioctahedral), calcite and zeolites (analcime and scolecite). In summary, Ca-Mg-Fe carbonates formed in stage 2, whereas calcite (±ankerite), zeolites and increasing amounts of smectites formed in stages 3 and 4. The main consequence of increasing the initial pCO2 in the simulations, was to increase the amount of dissolved basaltic glass to reach a certain pH value of one of the four different stages. A comparison of the simulated pH for the system initially saturated at a pCO2 of 20 bar vs. the one at 100 bar (Fig. 1), indicates that the former reached stage 3 after dissolution of ~100 g of basaltic glass, whereas the latter needed >200 g of dissolved basaltic glass to reach these pH values. This is in line with the lower concentration of protons that need to be neutralized by the basaltic glass at lower pCO2 (i.e. due to a decrease in proton activity), which results in corresponding higher final pH values for a given amount of dissolved glass. Comparison of the mineralogy between these two simulations (Fig. 2), indicates that at lower pCO2 (i.e. 20 bar) the CO2 mineralized dominantly into calcite, whereas other divalent cations (Mg and Fe) precipitated into di- and trioctahedral smectites. At elevated pCO2 (i.e. 100 bar), the CO2 mineralized dominantly into Ca-Mg-Fe carbonates with only minor amounts of smectites that precipitated.
The stages mentioned above are also clearly defined by the simulated aqueous CO2 species, proton activity and the total dissolved cation concentrations (Figs. 3 and 4 ). At low initial pCO2 (i.e. 20 bar), CO2(aq) was the major species in stages 1 and 2, and ionized to HCO3− with increasing pH (i.e. decreasing proton activity). The CO2 speciation is also consistent with the highest mobility of Ca, Mg and Fe defined by an increase in their total dissolved concentrations. Once the pH reached values of ~5–5.5 (stage 3), HCO3− became the main aqueous CO2 species followed by a rapid decrease in Ca, Mg, Fe and K concentrations associated with the precipitation of smectites. An increase in total dissolved Al, indicates that in stages 3 and 4, the mobility of Al influenced considerably the solubility of the smectites and kaolinite. The availability of Al increased with pH due to the increased stability of the Al(OH)4− species , which affected the solubility of kaolinite. Comparison of these observations with the simulations at elevated pCO2 (i.e. 100 bar), indicates that adding CO2 to the system increases the buffer capacity of the acid (i.e. CO2), which increases the stability of CO2(aq) and the total dissolved Ca, Mg and Fe concentrations. This increased acid buffer capacity with increasing pCO2 is also consistent with the findings of , who calculated the buffer capacity of CO2-rich solutions associated to mineral carbonation. In addition, the simulations reveal that addition of CO2 will enhance the mobility of Ca, Mg and Fe for the mineralization into Ca–Mg–Fe carbonates (i.e. Mg-Fe carbonate solid solutions, dolomite and ankerite) during stage 2, whereas reducing CO2 extends stages 3 and 4 where Ca, Mg and Fe become less mobile due to the availability of Al to precipitate smectites. The titration model indicates that the maximum amount of mineralized CO2, which corresponds to the steepest slope in a diagram of moles CO2 mineralized vs. grams of basalt added (Fig. 5), occurs in stage 2 of the reaction path.
The previous simulations did not take into account the dissolution kinetics of the basaltic glass, which can be used to determine the relative timeframe of the reaction paths described above. In the following, I will focus on the basaltic glass dissolution kinetics (i.e. omitting secondary mineral precipitation/dissolution kinetics) for the sake of delineating differences with the titration model. Since the dissolution rate of basaltic glass (according to Equations 5–7) is affected by pH and Al concentrations , a series of simulations have been carried out to test these effects on the formation of secondary alteration minerals. The simulated pH values (Fig. 6) show that after 500 days the systems saturated with an initial pCO2 of 40 and 100 bar had a pH buffered at ~5.5–6, which corresponds to stage 2 of the titration model (Fig. 1), whereas the system with an initial pCO2 of 20 bar, had a pH buffered at ~8–8.5, which corresponds to stage 3 of the titration model (Fig. 1). The overall sequence of simulated mineral assemblages of the kinetic model (Fig. 7) is similar to the titration model (Fig. 2), where stage 2 minerals consist mainly of chalcedony, kaolinite, and Ca-Mg-Fe carbonates and stage 3 of chalcedony, kaolinite, and Ca-(Fe) carbonates and smectites. The kinetic simulations indicate a rapid mineralization of CO2 and the dissolution of all the basaltic glass shortly after 200 days in both simulations (Figs. 7a and b). Therefore, modifying the initial pCO2 had a major effect on the secondary minerals formed but only a minor effect on the total amount of dissolved basaltic glass. Equation 6 indicates that the basaltic glass forward dissolution rate is very sensitive to varying the reactive surface area of the glass. For the simulations at pCO2 of 20 bar, varying the initial geometric surface area from 250 cm2/g to 50 cm2/g had a drastic effect on the simulated basaltic glass dissolution and secondary minerals formed (Fig. 7b and c). Reducing the reactive surface area was capable of considerably shifting the simulated mineral assemblages from stage 3 to stage 2 with the formation of Ca–Mg–Fe carbonates instead of smectites. In addition, the amount of dissolved basaltic glass per unit time was reduced, which led to a decrease in the total amount of carbonates formed. These simulations show the importance of evaluating the reactive surface area upon fluid-rock reaction, and may also point to the need to consider effects of mineral coatings on the primary mineral surfaces, and their potential to decrease mineral dissolution rates , . Finally, the simulated changes of rock volumes, and the contribution of basaltic glass, chalcedony, kaolinite, smectites and carbonates to the total volume is shown in Fig. 8. The results of these simulations indicate that the total volume of precipitated solids was higher in the model with an initial pCO2 of 100 bar, caused by the mineralization of higher amounts of chalcedony and kaolinite in comparison to the 20 bar pCO2 model, where the amount of smectite mineralization was more significant. These new findings indicate that precipitation of kaolinite and chalcedony could potentially lead to a clogging of pore space in a basaltic rock formation depending on the initial system porosity.
The titration and kinetic models were useful for determining the controlling factors for the mineralization of carbonates vs. smectites, and for defining ways to optimize the geochemical conditions for CO2 sequestration in basaltic rock formations. I will emphasize here some of the limitations of these models that may affect chemical mass transfers between water and rock in these simulations. One limitation was to assume that all the minerals are at partial equilibrium in this thermodynamic system, except for the dissolution rate of basaltic glass. At low temperatures (<100°C), mineral solubilities are affected by reaction kinetics, such as for example the precipitation of certain carbonates (i.e. magnesite, ankerite and dolomite) and smectites , , . Indeed, low temperature laboratory CO2-water-basalt experiments have shown that many of the secondary minerals formed at 40°C were poorly crystalline to amorphous phases . Currently, most kinetic models rely on transition state theory , which may not predict adequately measured laboratory mineral growth rates. In addition, it is difficult to reconcile measured laboratory mineral dissolution/precipitation rates with field observations . Despite these difficulties, laboratory CO2-water-basalt experiments have demonstrated that chalcedony, Ca–Fe- and Ca–Mg-carbonates mineralized within 124 days at 75°C and a pCO2 of ~25 bar . This is in line with the alteration mineralogy (i.e. quartz, Fe–Mg carbonate solid solutions and dolomite) observed in the Disko-Nuussuaq Island (West Greenland), resulting from the interaction of a CO2-rich water with basaltic rock formation at estimated temperatures of 80–126°C . These observations indicate that the simplified titration and kinetic models yield reasonable results with respect to the simulated alteration mineralogy. Caution should be taken with minerals like dolomite, which are predicted to precipitate in the simulations, but could occur as Ca–Mg–Fe carbonate solid solutions in the real system with a different structure and composition than the dolomite phase. Nevertheless, the present study did take into account the formation of Mg–Fe carbonate solid solutions, which were predicted to precipitate instead of magnesite and siderite, and were also observed to form in several different types of low temperature laboratory experiments , , , . Future simulations should also take into account that ankerite does not occur in nature, which may indicate the formation of a ferroan calcite instead.
Comparison of GEM vs. LMA methods
Simulating the interaction of CO2-water-basalt using the GEMS code package resulted in several improvements over our previous study  using PHREEQC based on the LMA method . The decisive advantages of the GEM method include: (i) the thermodynamic database structure, (ii) the possibility to include multicomponent and multisite mineral solid solutions, (iii) the definition of equations of state for gases and their mixtures (i.e. to calculate fugacity coefficients for non-ideal gases). The GEMS code package permits to enter all the standard thermodynamic properties (i.e. dG, dH, S, Cp and V functions) for all minerals, aqueous species and gases, including equation of state parameters. This enables to correct the thermodynamic properties for all species as a function of pressure and temperature , . Updates of the SUPCRT92 dataset  were successfully integrated in the MINES thermodynamic database with more recent thermodynamic data for aqueous species based on experimental work, and the mineral properties from the internally consistent Holland and Powell  database (Supplementary Materials). This database structure avoids the necessity to recalculate all the equilibrium reactions for mineral dissolution reactions as would be required by databases compiled for LMA type codes . The latter also require to define a set of master species (e.g. Al3+), which are used to calculate all subsequent equilibrium constant for every mineral dissolution reaction , . This may present certain disadvantages for simulating a reaction path covering a large range of pH values, where for example, the calculated activities of these master species may become insignificantly low. This may consequently lead to larger uncertainties in calculated reaction quotients for mineral dissolution reactions, which will strongly depend on the quality of selected thermodynamic data.
Comparison of the overall simulated pH evolution and mineralogical changes, in both, GEM and previous LMA study (Table 2), indicate consistent results but also fine differences. The new simulations using the GEM method yields a more detailed pH range and associated secondary mineral assemblage with reaction progress in comparison to the previous  simulations. The major change in the present study, was the capability to simulate the stability of tri-/dioctahedral smectite solid solution with mixing on two crystallographic sites, i.e. eight end members for dioctahedral smectites and five end members for trioctahedral smectites (Supplementary Materials Table 1). In addition, the simulations considered Ca–Mg–Fe carbonate solid solutions, without having a priori fixed compositions as was required by the LMA code. This capability can be exploited by the GEM method, because the chemical potentials of each component and phases are calculated to minimize the Gibbs energy of the system , . This also permits to calculate pH and the redox state of the system using the Gibbs energy minimization . In LMA codes, the pH is either entered or calculated based on charge balance, and the redox state of the system must be entered based on either a selected redox pairs [e.g. Fe(II)/Fe(III)], or a given initial Eh value. One of the drawbacks of the GEM method, was the slower computing times compared to LMA codes , which must be considered for simulating a complex interplay of mineral dissolution and precipitation kinetics . Nevertheless, recent efforts have been made to combine the advantages of both, GEM and LMA methods, leading to improved computing times for solving complex multicomponent and multphase systems , .
|This study (GEM method)||Gysi and Stefánsson |
|pH||~5–5.5||~4.5–6 (CO2-basalt alteration buffered)|
|Mineralogy||chalc, kln, Ca–Mg–Fe carb||chalc, Ca–(Mg)–Fe(II)–Fe(III) sme2, Al–Si minerals, Ca–Mg–Fe carb|
|pH||~8–8.5||>8 (basalt alteration buffered)|
|Mineralogy||chalc, kln, di/trioctahedral sme1, Ca–(Fe)-carb||Ca–Mg–Fe(II) sme2, Fe oxyhydroxides, zeo, cc|
|Mineralogy||chalc, kln, di/trioctahedral sme1, cc, zeo|
chalc, chalcedony; carb, carbonates; kln, kaolinite; sme, smectites; cc, calcite; zeo, zeolites.
1The composition of di- and trioctahedral smectites was allowed to vary in the simulations assuming an ideal solid solution model and mixing on two different crystallographic sites (Supplementary Materials Table 1).
2The composition of di- and trioctahedral smectites was fixed at the beginning of the simulations, and their solubilities were calculated according to their saturation in solution.
Challenges for optimizing CO2 sequestration in basaltic rock formations
The overall simulated reaction paths yield an idea on how to optimize the sequestration of CO2. Monitoring the borehole geochemistry during injection of CO2 into basaltic rock formations could be used together with Figs. 1–4 to determine the expected mineral assemblages forming under these conditions in the subsurface. A combination of pH, major CO2 aqueous species, and total dissolved element concentrations can be used as a geochemical signal for evaluating the reaction progress (stages 1–4) of the overall CO2-water-basalt interaction process. The amount of injected CO2 can be fine tuned to achieve a higher mobility of Ca, Mg and Fe for carbonate mineralization, which is favorable in stage 2, and to decrease the mobility of Al to avoid smectite formation. To avoid clogging of pore spaces (Fig. 8), however, the amount of CO2 injected should be kept below a certain threshold, depending on the initial basalt porosity and the flow pathways of the fluid. The total simulated volume increased by up to ~30% of the total initial rock volume for the simulations at pCO2 of 100 bar, whereas it increased only by up to ~20% for the simulations at pCO2 of 20 bar in the kinetic simulations (Fig. 8). This highlights some of the limitations of the simple titration and kinetic models, and the importance of studying permeability and porosity changes associated to fluid-rock reactions for delineating how to maximize the mineralization of CO2 using reactive mass transport simulations. One of the future challenges to address, will include a correct estimation of volume changes associated to the formation of different types of phyllosilicates, including hydration of the interlayer that will lead to larger occupied pore space volumes , , and consequently reduce the reactive surface area of the basaltic rock. Dilution of the injected CO2-rich waters with aquifer groundwaters, and also variations in the composition of the injected fluids may lead to more complexity in the understanding of the fluid-rock interaction process. To interpret correctly field data, will require a combination of hydrological tracers and reactive transport simulations to determine various mixing models . In addition, the effects of using solid solutions with fixed composition (LMA codes) vs. variable compositions (GEM codes) need to be carefully investigated, to delineate the effects on the overall calculated reaction paths and mass balances. It has been shown that smectites formed through weathering of crystalline basalts display a composition similar to amorphous phases formed by the alteration of basaltic glass on geological timescales . The basalt crystallinity may have significant effects on the amounts of Ca, Mg and Fe available for carbonate formation during CO2 injection into a basaltic aquifer on shorter timescales. Recent laboratory experiments have shown that the release rate of Mg and Fe over Ca is different for crystalline basalts (i.e. plagioclase, pyroxene, and olivine) vs. basaltic glass at low pH , and future numerical modeling studies should consider these kinetic effects to test and apply the simulations to different types of basaltic rock formations.
Conclusions and outlook
In this study, I have shown that using the GEM method for simulating water-rock interaction permits to model in finer details chemical reactions, stable mineralogy and pH for the CO2-water-basalt reaction paths. Addition of multicomponent and multisite mineral solid solutions for phyllosilicates, zeolites and carbonates is an important step forward for simulating this type of system. Difficulties in retrieving experimental data for phyllosilicates, indicates a need to improve and test theoretical models (i.e. oxide summation technique) for estimating their thermodynamic properties. Further, the simulated titration and kinetic models may be used to guide future reactive mass transport simulations, that will permit to link the chemical model of the present study with physical models that considers porosity, permability and reactive surface area changes upon CO2-water-basalt interaction in basaltic rock formations. The latter can be achieved by coupling GEMS with reactive transport codes , . Finally, the GEM and LMA methods have recently been combined into a new extended xLMA model , which may together with future updates of the MINES database, considerably improve our future capabilities to simulate fluid-rock reactions in the subsurface.
A collection of invited papers based on presentations at the International Symposium on Solubility Phenomena and Related Equilibrium Processes (ISSP-17), Geneva, 24–29 July 2016.
This study was supported by the Colorado School of Mines with a professional development grant to A. P. Gysi. This manuscript benefited from discussions with D. Kulik for the implementation of the kinetics and mineral solid solution models in GEMS, and was significantly improved by the thoughtful reviews of an anonymous reviewer. I would like to thank editor H. Burrows for organizing this special issue, and the organizers of the ISSP17 conference, M. Filella and W. Hummel, for inviting me to give a lecture in Geneva, Switzerland.
 E. H. Oelkers, S. R. Gislason, J. Matter. Elements. 4, 333 (2008).10.2113/gselements.4.5.333Search in Google Scholar
 S. R. Gislason, E. H. Oelkers. Science344, 373 (2014).10.1126/science.1250828Search in Google Scholar PubMed
 J. M. Matter, P. B. Kelemen. Nat. Geosci.2, 837 (2009).10.1038/ngeo683Search in Google Scholar
 B. P. McGrail, H. T. Schaef, A. M. Ho, Y.-J. Chien, J. J. Dooley, C. L. Davidson. J. Geophys. Res.111, 1 (2006).Search in Google Scholar
 J. M. Matter, W. S. Broecker, M. Stute, S. R. Gislason, E. H. Oelkers, A. Stefánsson, D. Wolff-Boenisch, E. Gunnlaugsson, G. Axelsson, G. Björnsson. Energy Procedia. 1, 3641 (2009).10.1016/j.egypro.2009.02.160Search in Google Scholar
 J. M. Matter, M. Stute, S. Ó. Snaebjörnsdottir, E. H. Oelkers, S. R. Gislason, E. S. Aradottir, B. Sigfusson, I. Gunnarsson, H. Sigurdardottir, E. Gunnlaugsson, G. Axelsson, H. A. Alfredsson, D. Wolff-Boenisch, K. Mesfin, D. Fernandez de la Reguera Taya, J. Hall, K. Dideriksen, W. S. Broecker. Science352, 1312 (2016).10.1126/science.aad8132Search in Google Scholar PubMed
 P. S. Neuhoff, T. Fridriksson, D. K. Bird. Am. J. Sci.299, 467 (1999).10.2475/ajs.299.6.467Search in Google Scholar
 K. L. Rogers, P. S. Neuhoff, A. K. Pedersen, D. K. Bird. Lithos. 92, 55 (2006).10.1016/j.lithos.2006.04.002Search in Google Scholar
 P. S. Neuhoff, K. L. Rogers, L. S. Stannius, D. K. Bird, A. K. Pedersen. Lithos. 92, 33 (2006).10.1016/j.lithos.2006.03.028Search in Google Scholar
 A. P. Gysi, A. Stefánsson. Geochim. Cosmochim. Acta. 81, 129 (2012).10.1016/j.gca.2011.12.012Search in Google Scholar
 A. P. Gysi, A. Stefánsson. Chem. Geol.306–307, 146 (2012).10.1016/j.chemgeo.2012.03.006Search in Google Scholar
 T. H. Van Pham, P. Aagaard, H. Hellevang. Geochem. Trans.13, 5 (2012).10.1186/1467-4866-13-5Search in Google Scholar
 A. P. Gysi, A. Stefánsson. Geochim. Cosmochim. Acta75, 4728 (2011).10.1016/j.gca.2011.05.037Search in Google Scholar
 E. S. P. Aradóttir, E. L. Sonnenthal, G. Björnsson, H. Jónsson. Int. J. Greenh. Gas Control. 9, 24 (2012).10.1016/j.ijggc.2012.02.006Search in Google Scholar
 D. L. Parkhurst, C. A. J. Appelo. U. S. Geol. Surv. Tech. Meth.A43, 497 (2013).Search in Google Scholar
 T. Xu, N. Spycher, E. Sonnenthal, G. Zhang, L. Zheng, K. Pruess. Comput. Geosci.37, 763 (2011).10.1016/j.cageo.2010.10.007Search in Google Scholar
 E. S. P. Aradóttir, E. L. Sonnenthal, H. Jónsson. Chem. Geol.304–305, 26 (2012).10.1016/j.chemgeo.2012.01.031Search in Google Scholar
 K. M. Krupka, K. J. Cantrell, B. P. McGrail. U.S. Department of Energy report PNNL-19766 (2010).Search in Google Scholar
 C. Steefel, D. DePaolo, P. Lichtner. Earth Planet. Sci. Lett.240, 539 (2005).10.1016/j.epsl.2005.09.017Search in Google Scholar
 H. Helgeson, R. Garrels, F. MacKenzie. Geochim. Cosmochim. Acta. 33, 455 (1969).10.1016/0016-7037(69)90127-6Search in Google Scholar
 C. Zhu, P. Lu, Z. Zheng, J. Ganor. Geochim. Cosmochim. Acta. 74, 3963 (2010).10.1016/j.gca.2010.04.012Search in Google Scholar
 R. M. Garrels. Genesis of Some Ground Waters from Igneous Rocks, Wiley, New York (1967).Search in Google Scholar
 H. Helgeson. Geochim. Cosmochim. Acta. 32, 853 (1968).10.1016/0016-7037(68)90100-2Search in Google Scholar
 D. Wolff-Boenisch, S. R. Gislason, E. H. Oelkers, C. V. Putnis. Geochim. Cosmochim. Acta. 68, 4843 (2004).10.1016/j.gca.2004.05.027Search in Google Scholar
 S. R. Gislason, E. H. Oelkers. Geochim. Cosmochim. Acta. 67, 3817 (2003).10.1016/S0016-7037(03)00176-5Search in Google Scholar
 L. Marini. Geological Sequestration of Carbon Dioxide: Thermodynamics, Kinetics, and Reaction Path Modeling, Elsevier, Amsterdam (2007).Search in Google Scholar
 V. Daux, C. Guy, T. Advocat, J.-L. Crovisier, P. Stille. Chem. Geol.142, 109 (1997).10.1016/S0009-2541(97)00079-XSearch in Google Scholar
 G. Berger, C. Claparols, C. Guy, V. Daux. Geochim. Cosmochim. Acta. 58, 4875 (1994).10.1016/0016-7037(94)90218-6Search in Google Scholar
 P. Aagaard, H. C. Helgeson. Am. J. Sci.282, 237 (1982).10.2475/ajs.282.3.237Search in Google Scholar
 D. A. Kulik, T. Wagner, S. V. Dmytrieva, G. Kosakowski, F. F. Hingerl, K. V. Chudnenko, U. R. Berner. Comput. Geosci.17, 1 (2013).Search in Google Scholar
 A. P. Gysi, A. Stefánsson. Chem. Geol.306–307, 10 (2012).10.1016/j.chemgeo.2012.02.016Search in Google Scholar
 G. J. Stockmann, D. Wolff-Boenisch, S. R. Gislason, E. H. Oelkers. Chem. Geol.284, 306 (2011).10.1016/j.chemgeo.2011.03.010Search in Google Scholar
 A. P. Gysi, A. Stefánsson. Mineral. Mag.72, 55 (2008).10.1180/minmag.2008.072.1.55Search in Google Scholar
 T. Holland, R. Powell. J. Metamorph. Geol.16, 309 (1998).10.1111/j.1525-1314.1998.00140.xSearch in Google Scholar
 R. A. Robie, B. S. Hemingway. U. S. Geol. Surv. Bull.2131, 458 (1995).Search in Google Scholar
 D. A. Sverjensky, E. L. Shock, H. C. Helgeson. Geochim. Cosmochim. Acta. 61, 1359 (1997).10.1016/S0016-7037(97)00009-4Search in Google Scholar
 J. C. Tanger, H. C. Helgeson. Am. J. Sci. 288, 19 (1988).10.2475/ajs.288.1.19Search in Google Scholar
 H. C. Helgeson, D. H. Kirkham, G. C. Flowers. Am. J. Sci.281, 1249 (1981).10.2475/ajs.281.10.1249Search in Google Scholar
 E. L. Shock, H. C. Helgeson, D. A. Sverjensky. Geochim. Cosmochim. Acta. 53, 2157 (1989).10.1016/0016-7037(89)90341-4Search in Google Scholar
 J. W. Johnson, E. H. Oelkers, H. C. Helgeson. Comput. Geosci.18, 899 (1992).10.1016/0098-3004(92)90029-QSearch in Google Scholar
 G. D. Miron, T. Wagner, D. A. Kulik, C. A. Heinrich. Geochim. Cosmochim. Acta. 187, 41 (2016).10.1016/j.gca.2016.04.026Search in Google Scholar
 G. D. Miron, D. A. Kulik, S. V. Dmytrieva, T. Wagner. Appl. Geochem.55, 28 (2015).10.1016/j.apgeochem.2014.10.013Search in Google Scholar
 B. Tagirov, J. Schott. Geochim. Cosmochim. Acta. 65, 3965 (2001).10.1016/S0016-7037(01)00705-0Search in Google Scholar
 D. A. Palmer, P. Bénézeth, D. J. Wesolowski. Geochim. Cosmochim. Acta. 65, 2081 (2001).10.1016/S0016-7037(01)00584-1Search in Google Scholar
 P. Bénézeth, D. A. Palmer, D. J. Wesolowski. Geochim. Cosmochim. Acta. 65, 2097 (2001).10.1016/S0016-7037(01)00585-3Search in Google Scholar
 P. S. Neuhoff, thesis, Stanford University (2000).Search in Google Scholar
 T. Wagner, D. A. Kulik, F. F. Hingerl, S. V. Dmytrieva. Can. Mineral.50, 1173 (2012).10.3749/canmin.50.5.1173Search in Google Scholar
 R. A. Robinson, R. H. Stokes. Electrolyte solutions, Butterworths, London (1968).Search in Google Scholar
 H. C. Helgeson, D. H. Kirkham. Am. J. Sci.274, 1199 (1974).10.2475/ajs.274.10.1199Search in Google Scholar
 R. Stryjek, J. H. Vera. Can. J. Chem. Eng.64, 323 (1986).10.1002/cjce.5450640224Search in Google Scholar
 P. Proust, J. H. Vera. Can. J. Chem. Eng.67, 170 (1989).10.1002/cjce.5450670125Search in Google Scholar
 J. A. Chermak, J. D. Rimstidt. Am. Mineral.74, 1023 (1989).Search in Google Scholar
 T. J. B. Holland, Am. Mineral.74, 5 (1989).Search in Google Scholar
 R. Berman, T. Brown. Contrib. Mineral. Petrol.89, 168 (1985).10.1007/BF00379451Search in Google Scholar
 D. Wolff-Boenisch. Energy Procedia. 4, 3738 (2011).10.1016/j.egypro.2011.02.307Search in Google Scholar
 G. D. Saldi, D. Daval, G. Morvan, K. G. Knauss. Geochim. Cosmochim. Acta. 118, 157 (2013).10.1016/j.gca.2013.04.029Search in Google Scholar
 G. D. Saldi, J. Schott, O. S. Pokrovsky, Q. Gautier, E. H. Oelkers. Geochim. Cosmochim. Acta. 83, 93 (2012).10.1016/j.gca.2011.12.005Search in Google Scholar
 O. S. Pokrovsky, S. V. Golubev, J. Schott, A. Castillo. Chem. Geol.260, 317 (2009).Search in Google Scholar
 H. Hellevang, V. T. H. Pham, P. Aagaard. Int. J. Greenh. Gas Control. 15, 3 (2013).10.1016/j.ijggc.2013.01.027Search in Google Scholar
 A. F. White, S. L. Brantley. Chem. Geol.202, 479 (2003).10.1016/j.chemgeo.2003.03.001Search in Google Scholar
 C. S. Romanek, C. Jiménez-López, A. R. Navarro, M. Sanchez-Román, N. Sahai, M. Coleman. Geochim. Cosmochim. Acta. 73, 5361 (2009).10.1016/j.gca.2009.05.065Search in Google Scholar
 C. M. Bethke. Geochemical and Biogeochemical Reaction Modeling, Cambridge University Press, New York (2007).10.1017/CBO9780511619670Search in Google Scholar
 I. K. Karpov, K. V. Chudnenko, D. A. Kulik, O. V. Avchenko, V. A. Bychinskii. Geochemistry Int.39, 1108 (2001).Search in Google Scholar
 A. M. M. Leal, D. A. Kulik, G. Kosakowski. Adv. Water Resour.88, 231 (2016).10.1016/j.advwatres.2015.11.021Search in Google Scholar
 A. M. M. Leal, D. A. Kulik, G. Kosakowski, M. O. Saar. Adv. Water Resour.96, 405 (2016).10.1016/j.advwatres.2016.08.008Search in Google Scholar
 O. Vidal, B. Dubacq. Geochim. Cosmochim. Acta731, 6544 (2009).10.1016/j.gca.2009.07.035Search in Google Scholar
 B. M. J. Thien, G. Kosakowski, D. A. Kulik. Geotherm. Energy. 3, 11 (2015).10.1186/s40517-015-0031-7Search in Google Scholar
 J. L. Crovisier, V. Daux. Chem. Geol.84, 261 (1990).10.1016/0009-2541(90)90232-VSearch in Google Scholar
 S. Gudbrandsson, D. Wolff-Boenisch, S. R. Gislason, E. H. Oelkers. Geochim. Cosmochim. Acta. 75, 5496 (2011).10.1016/j.gca.2011.06.035Search in Google Scholar
The online version of this article (DOI: 10.1515/pac-2016-1016) offers supplementary material, available to authorized users.
©2017 IUPAC & De Gruyter. This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License. For more information, please visit: http://creativecommons.org/licenses/by-nc-nd/4.0/