Numerical simulations of CO₂ sequestration in basaltic rock formations: challenges for optimizing mineral-fluid reactions

Reaction path model

Gaseous CO₂ dissolves in natural waters to form carbonic acid, which ionizes to bicarbonate at mildly acidic pH conditions according to,

where H₂CO₃* represents the species CO₂(aq) and minor H₂CO₃ (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. Ca²⁺, Mg²⁺ and Fe²⁺) for the mineralization of CO₂ 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 Ca²⁺, and the phyllosilicates that may compete for Ca²⁺, Mg²⁺ and Fe²⁺. 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 [20], [21]. This overall rock dissolution process leads to a pH increase similar to an acid-base titration [22], where the basalt will act as the base by consuming protons and releasing cations upon dissolution, and the water with its acids such as CO₂ 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 [20], [23].

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 [24], [25], and are therefore ideal for prospective CO₂ 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 [13], [26], 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 [25],

where A_A is a pre-exponential factor including the rate constant (10^−5.6 mol of Si/cm²/s), E_A an activation energy of 25.5 kJ/mol, R the ideal gas constant, T temperature in K, and a_i the activity of the ith aqueous species. The overall forward dissolution rate (r₊) based on the geometric reactive surface area (s_geo,t in cm²) 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 s_geo is the initial reactive surface area in cm², m_t the mole amount at time t and m₀ 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 Ca²⁺ and Na⁺, whereas higher charged cations such as Si⁴⁺ and Al³⁺ will form stronger bonds within the glass structure [25]. 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 [24], [25], [27], [28]. 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 [29], where Q is the reaction quotient and K the equilibrium constant for the solubility of the basaltic glass.

Reaction path modeling

Numerical simulations of the CO₂-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 [30]. These simulations are different than in our previous study [13], 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 [24], [25], [31], [10], [32]. The simplified chemical formula of this natural basaltic glass can be described by: K_0.008Na_0.08Ca_0.27Mg_0.26Fe^(II)_0.169Fe^(III)_0.012SiAl_0.35O_3.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 [31], [10], [11] and the LMA simulations of [13], [33]. 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 [13]. Its initial composition was varied by saturating it at 50°C with different initial CO₂ gas partial pressures (pCO₂) 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 CO₂-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 CO₂-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 cm²/g [25], which corresponds to spherical particles with a radius of 45–125 μm. In both of these models, the effects of changing pCO₂ 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).

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 [30]. 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 [34], [35], and aqueous species from SUPCRT92 (slop98.dat; [36], [37], [38], [39], [40]), and was updated with the internally consistent dataset of [41] for K, Na and Si aqueous species optimized using the GEMSFITS code package [42], and with more recent experimental data for Al aqueous species [43], [44], [45]. The properties of zeolites were derived from the study of [46], and the properties of phyllosilicates (i.e. saponite, beidellite and nontronite) using the approach from the study of [13] as described below. The standard state adopted in this study was unit activity for pure H₂O 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 [47]. The activity model used for aqueous species was the extended Debye-Hückel according to [48] and using the Debye-Hückel parameters of [38], [49]. The equation of state used for real gases (i.e. CO₂–CH₄–CO–H₂O) was the Peng-Robinson-Stryjek-Vera (PRSV) model [50], [51].

The thermodynamic properties of di- and trioctahedral smectite end members were estimated using the combined approach of [13], 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 [52], the entropy according to [53], and the heat capacity function according to [54]. 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 [13], 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 [52], [53], [54].

Titration model

Results of the simulations of basaltic glass reacted at 50°C with the CO₂-rich aqueous solutions saturated initially at pCO₂ 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 CO₂ (aq) leading to the formation of bicarbonate (HCO₃⁻) 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 (SiO₂), 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 (SiO₂), 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 pCO₂ 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 pCO₂ of 20 bar. At these conditions, the stable mineral assemblage consisted of chalcedony (SiO₂), 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 pCO₂ 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 pCO₂ 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 pCO₂ (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 pCO₂ (i.e. 20 bar) the CO₂ mineralized dominantly into calcite, whereas other divalent cations (Mg and Fe) precipitated into di- and trioctahedral smectites. At elevated pCO₂ (i.e. 100 bar), the CO₂ mineralized dominantly into Ca-Mg-Fe carbonates with only minor amounts of smectites that precipitated.

Fig. 1:

Titration model showing the simulated pH as a function of basaltic glass added at 50°C and pCO₂ of 20, 40 and 100 bar. The numbers indicate different pH buffered conditions: (1) CO₂-water-buffered; (2) SiO₂+kaolinite+Ca-Mg-Fe carbonate solid solutions; (3) SiO₂+kaolinite+smectites+Ca-(Fe)-carbonates; (4) SiO₂+kaolinite+smectites+zeolites+calcite.

Fig. 2:

Titration model showing the simulated secondary mineralogy as a function of basaltic glass added at 50°C and pCO₂ of (a) 20, (b) 40 and (c) 100 bar. The numbers indicate different pH buffered conditions (see also Fig. 1): (1) CO₂-water-buffered; (2) SiO₂+kaolinite+Ca-Mg-Fe carbonate solid solutions; (3) SiO₂+kaolinite+smectites+Ca-(Fe)-carbonates; (4) SiO₂+kaolinite+smectites+zeolites+calcite.

The stages mentioned above are also clearly defined by the simulated aqueous CO₂ species, proton activity and the total dissolved cation concentrations (Figs. 3 and 4 ). At low initial pCO₂ (i.e. 20 bar), CO₂(aq) was the major species in stages 1 and 2, and ionized to HCO₃⁻ with increasing pH (i.e. decreasing proton activity). The CO₂ 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), HCO₃⁻ became the main aqueous CO₂ 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)₄⁻ species [45], which affected the solubility of kaolinite. Comparison of these observations with the simulations at elevated pCO₂ (i.e. 100 bar), indicates that adding CO₂ to the system increases the buffer capacity of the acid (i.e. CO₂), which increases the stability of CO₂(aq) and the total dissolved Ca, Mg and Fe concentrations. This increased acid buffer capacity with increasing pCO₂ is also consistent with the findings of [55], who calculated the buffer capacity of CO₂-rich solutions associated to mineral carbonation. In addition, the simulations reveal that addition of CO₂ 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 CO₂ 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 CO₂, which corresponds to the steepest slope in a diagram of moles CO₂ mineralized vs. grams of basalt added (Fig. 5), occurs in stage 2 of the reaction path.

Fig. 3:

Titration model showing the simulated aqueous speciation and total dissolved element concentrations as a function of basaltic glass added at 50°C and pCO₂ of 20 bar. The numbers indicate different stages corresponding to mineral assemblages buffering the pH and controlling the element mobility (Figs. 1 and 2). Stage 3 (SiO₂+kaolinite+smectites+Ca-(Fe)-carbonates) is the most extensive and marks the appearance of smectites and a decrease in Ca, Mg and Fe mobility.

Fig. 4:

Titration model showing the simulated aqueous speciation and total dissolved element concentrations as a function of basaltic glass added at 50°C and pCO₂ of 100 bar. The numbers indicate different stages corresponding to mineral assemblages buffering the pH and controlling the element mobility (Figs. 1 and 2). The increased CO₂ concentrations in these simulations lead to an increase of the mobility of Ca, Mg, and Fe, and stage 2 (SiO₂+kaolinite+Ca-Mg-Fe carbonate solid solutions) becoming more extensive in comparison to the lower CO₂ system (Fig. 3).

Fig. 5:

Titration model showing the total amounts of mineralized CO₂ as a function of basaltic glass added at 50°C and pCO₂ of 20, 40, 100 bar. The numbers indicate the extent and controls of stage 2 (SiO₂+kaolinite+Ca–Mg–Fe carbonate solid solutions) on carbonate mineralization.

Kinetic model

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 [25], 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 pCO₂ 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 pCO₂ 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 CO₂ and the dissolution of all the basaltic glass shortly after 200 days in both simulations (Figs. 7a and b). Therefore, modifying the initial pCO₂ 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 pCO₂ of 20 bar, varying the initial geometric surface area from 250 cm²/g to 50 cm²/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 [32], [56]. 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 pCO₂ of 100 bar, caused by the mineralization of higher amounts of chalcedony and kaolinite in comparison to the 20 bar pCO₂ 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.

Fig. 6:

Kinetic model at 50°C showing the simulated pH as a function of reaction time (days) for the dissolution of 150 g basaltic glass and initial pCO₂ of 20, 40 and 100 bar.

Fig. 7:

Kinetic model at 50°C showing the simulated secondary mineralogy (mole) as a function of reaction time (days) for the dissolution of 150 g basaltic glass. In these simulations pCO₂ and the initial reactive geometric surface area (A_geo) were varied: (a) pCO₂ of 100 bar and A_geo of 250 cm²/g, (b) pCO₂ of 20 bar and A_geo of 250 cm²/g, and (c) pCO₂ of 20 bar and A_geo of 50 cm²/g.

Fig. 8:

Kinetic model at 50°C showing the simulated mineral volume changes (cm³) as a function of reaction time (days) for the dissolution of 150 g basaltic glass at (a) pCO₂ of 20 bar and (b) pCO₂ of 100 bar.

Model limitations

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 CO₂ 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 [57], [58], [59]. Indeed, low temperature laboratory CO₂-water-basalt experiments have shown that many of the secondary minerals formed at 40°C were poorly crystalline to amorphous phases [10]. Currently, most kinetic models rely on transition state theory [29], 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 [60]. Despite these difficulties, laboratory CO₂-water-basalt experiments have demonstrated that chalcedony, Ca–Fe- and Ca–Mg-carbonates mineralized within 124 days at 75°C and a pCO₂ of ~25 bar [11]. 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 CO₂-rich water with basaltic rock formation at estimated temperatures of 80–126°C [8]. 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 [10], [11], [56], [61]. 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 CO₂-water-basalt using the GEMS code package resulted in several improvements over our previous study [13] using PHREEQC based on the LMA method [15]. 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 [30], [47]. Updates of the SUPCRT92 dataset [40] 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 [34] 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 [13]. The latter also require to define a set of master species (e.g. Al³⁺), which are used to calculate all subsequent equilibrium constant for every mineral dissolution reaction [13], [62]. 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 [13] 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 [30], [63]. This also permits to calculate pH and the redox state of the system using the Gibbs energy minimization [63]. 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 [30], which must be considered for simulating a complex interplay of mineral dissolution and precipitation kinetics [59]. 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 [64], [65].

Challenges for optimizing CO₂ sequestration in basaltic rock formations

The overall simulated reaction paths yield an idea on how to optimize the sequestration of CO₂. Monitoring the borehole geochemistry during injection of CO₂ 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 CO₂ 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 CO₂-water-basalt interaction process. The amount of injected CO₂ 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 CO₂ 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 pCO₂ of 100 bar, whereas it increased only by up to ~20% for the simulations at pCO₂ 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 CO₂ 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 [66], [67], and consequently reduce the reactive surface area of the basaltic rock. Dilution of the injected CO₂-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 [14]. 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 [68]. The basalt crystallinity may have significant effects on the amounts of Ca, Mg and Fe available for carbonate formation during CO₂ 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 [69], and future numerical modeling studies should consider these kinetic effects to test and apply the simulations to different types of basaltic rock formations.