Unveiling the Essential Parameters Driving Mineral Reactions during CO₂ Storage in Carbonate Aquifers through Proxy Models

Featured Application

The findings of this study offer support for the design of laboratory experiments, monitoring indicators, and the selection of parameters for history matching, ultimately enhancing the accuracy of databases of mineral reactions induced by CO₂ injection.

Abstract

Numerical simulation is a commonly employed technique for studying carbon dioxide (CO₂) storage processes in porous media, particularly saline aquifers. It enables the representation of diverse trapping mechanisms and the assessment of CO₂ retention capacity within the subsurface. The intricate physicochemical phenomena involved necessitate the incorporation of multiphase flow, accurate depiction of fluid and rock properties, and their interactions. Among these factors, geochemical reaction rates and mechanisms are pivotal for successful CO₂ trapping in carbonate reactive rocks. However, research on kinetic parameters and the influence of lithology on CO₂ storage remains limited. This limitation is partly due to the challenges faced in laboratory experiments, where the time scale of the reactions and the lack of in situ conditions hinder accurate measurement of mineral reaction rates. This study employs proxy models constructed using response surfaces calibrated with simulation results to address uncertainties associated with geochemical reactions. Monte Carlo simulation is utilized to explore a broader range of parameters and identify influential factors affecting CO₂ mineralization. The findings indicate that an open database containing kinetic parameters can support uncertainty assessment. Additionally, the proxy models effectively represent objective functions related to CO₂ injectivity and mineralization, with calcite dissolution playing a predominant role. pH, calcite concentration, and CO₂ injection rate significantly impact dolomite precipitation, while quartz content remains unaffected.

1. Introduction

The use of fossil fuels as the dominant energy source in today’s energy mix has raised concerns about the impact of human-generated CO₂ emissions on global warming. If not addressed, these concerns could lead to irreversible consequences. Carbon Capture and Storage (CCS) has emerged as a potential solution to mitigate these concerns [1]. The ideal scenario for CO₂ sequestration would involve a large storage area and effective trapping mechanisms to prevent the migration of injected CO₂ back to the surface. This approach should also be cost effective and have minimal environmental impacts [2].

Numerical simulation is commonly used to study the CO₂ storage process in porous media as it can model various trapping mechanisms and estimate the capacity of subsurface CO₂ retention. The numerical simulation models can accurately capture the various complex physicochemical phenomena, including the effects of multiphase flow, rock–fluid interactions, CO₂ solubility and diffusion in water, relative permeabilities, capillary pressure hysteresis, mineralization, and acidic reactions, all of which are crucial for CO₂ trapping in porous media during CCS simulation.

The rates of geochemical reactions play a crucial role in the dissolution and mineralization of injected CO₂, which is essential to permanently trap the CO₂ and reduce the chances of leakage into the atmosphere. These reactions can occur shortly after injection and persist over time scales ranging from hundreds to thousands of years [3,4]. In order to ensure the long-term environmental sustainability of a CCS project, it is mandatory to have a comprehensive understanding of the rates and mechanisms of key geochemical reactions and their impacts on the performance of CO₂ storage and petrophysical properties of the reservoir.

The absence of in situ conditions in laboratory experiments, specifically the lack of pore structure contributions, as seen in real storage projects, also contributes to potential inaccuracies when measuring the kinetic parameters of mineral reactions. Overall, predicting in situ reaction rates remains highly complex [5,6,7,8,9] and has prompted researchers to use proxy [10,11] and batch models [12] or analytical approaches [13] to evaluate the geochemical impact.

In the realm of geochemical models coupled with flow models in porous media, the Transition State Theory (TST) [14] is frequently employed to simulate CO₂ mineralization [15,16,17,18,19,20,21,22]. However, the impact of kinetic parameters of TST on simulation outcomes has not been sufficiently researched, particularly when considering the potential influence of lithology on carbonate mineral dissolution and mineralization associated with CO₂ injection. Limited research has been conducted on the role of lithology in CO₂ storage, likely due to insufficient data related to mineral reactions for CO₂ sequestration. The data scarcity complicates the numerical assessment of sequestration and increases uncertainty [1,12]. The evaluation of the uncertainty involved is crucial for developing and implementing effective CCS strategies in geological formations that are not typically prioritized for carbon storage, such as carbonate rocks.

Therefore, this study aims to contribute to the understanding of carbonate reactions during CO₂ injection using robust reactive-transport models based on the reported geochemical data and on the TST geochemical model. Consequently, gaining a deeper understanding of the factors that influence carbonate reactions will facilitate the identification of appropriate storage targets based on rock mineralogy and key kinetic parameters. The proxy model and the workflow presented here will help optimize decision parameters (e.g., injection rate) based on uncertain parameters (e.g., rock mineralogy). Additionally, the conclusions of this study can provide valuable insight for designing laboratory experiments, helping to identify the essential properties to be analyzed and the most influential parameters to be varied during a history-matching of laboratory coreflood results with numerical models.

Therefore, the innovative aspects of this research can be summarized as follows:

• It goes beyond solely evaluating geological uncertainties, as seen in [23,24];
• It adopts a comprehensive approach to characterize the mechanisms of CO₂ trapping, encompassing solubility, residual (e.g., the trapped CO₂ saturation due to relative permeability and capillarity hysteresis), and mineral trapping. A particular focus is on understanding the reactions involving carbonate minerals and the crucial parameters that govern these reactions. This study uses a stochastic approach to model the kinetic parameters, which is different from other studies conducted in the literature, which either model the geochemical reactions deterministically [11,25,26] or neglect them altogether [27,28].

2. Methodology

This section will outline the geochemical and geological models employed in this study, describing each model and its corresponding characterization. Numerous simulation results (SR) obtained from those models generated by Latin Hypercube sampling are used to create computationally efficient proxy models based on quadratic response surfaces [29] calibrated with hundreds of numerical simulations, which work as training experiments. The proxy models are then used to generate thousands of new experiments by Monte Carlo simulations (MCS), thereby representing the SR for each objective function (OF) [30]. Therefore, MCS enables the exploration of a wider range of parameter values, and its results can be used to investigate the most influential parameters for CO₂ mineralization. Furthermore, the proxy models also enable a granular analysis of the impact of a parameter (and any modifications in its value) on the objective function (OF) [31].

The main assumptions In this paper are as follows:

• Geomechanical aspects and caprock behavior are not modeled;
• There is no occurrence of natural fractures in the carbonate rock matrix. This concern was addressed by Machado et al. in [22,32], concluding that calcite dissolution is the primary reaction. However, the effects of calcite dissolution can be partially offset by the precipitation of halite and dawsonite, particularly in the areas surrounding the wellbore and within natural fractures;
• Pure CO₂ stream is injected, e.g., without impurities and free water content;
• Dry-out effect due to water vaporization with CO₂ injection is not modeled;
• There is no change in wettability since the caprock, where wettability alterations may occur as a result of CO₂ adsorption [33,34], is not included in our model;
• This study focuses on the technical subsurface phenomena related to the uncertainty of CO₂ mineralization in porous media. It does not include the assessment of risks and uncertainties associated with surface processes [35], process efficiency [36,37], or economic aspects [38].

2.1. Modeling Geochemical Reactions

CO₂ injection for CCS was numerically simulated using CMG-GEM [39], version. 2023.30, considering the following features:

• Diffusion coefficient (D) for supercritical CO₂ in brine equals 3.65 × 10⁻⁵ cm²/s, according to Ahmadi et al. [40], with an average value over the pressure range considered during the simulations.
• CO₂ solubility in brine according to the method by Li and Nghiem [41] based on Henry’s law. This model is based on Henry’s constant calculation according to Equation (1) as a function of pressure and temperature. Still, the effect of salt on the gas solubility in the aqueous phase is modeled by the salting-out coefficient [42].

$$ln\left(H_i\right)=ln\left(H_i^{*}\right)+\frac{\overline{v_i}}{RT}\left(p-p^{*}\right)$$

(1)

where

$H_i$: Henry’s constant at current pressure (p) and temperature (T);

$H_i^{*}$: Henry’s constant at reference pressure (p^*) and temperature (T);

$\overline{v_i}$: partial molar volume at infinite dilution;

R: universal gas constant;

i: species dissolved in water (CO₂ in this work).

The$H_i^{*}$ constant was computed considering the initial conditions (prior to the commencement of CO₂ injection) of the models in

Table 1. Solubility parameter at the initial conditions of the reservoir model.

Properties	Values
Initial pressure	11,800 kPa @ 1200 m
Temperature	80 °C
Salinity	50,000 ppm (NaCl)
$H_i^{*}$	5.66 × 105

.

• Water acidic reactions for bicarbonate and carbonate ions generation, using kinetic parameters from the PHREEQC database [5,8]:

OH⁻_(aq) + H⁺_(aq) = H₂O_(aq)

(2)

CO_2(g) + H₂O_(aq) = H⁺_(aq) + HCO₃⁻_(aq)

(3)

CO₃²⁻_(aq) + H⁺_(aq) = HCO₃⁻_(aq)

(4)

The synthetic water composition (ions with non-zero initial concentration) assigned for this case is represented in

Table 2. Ionic composition of formation brine.

Ions	Concentration (mol/kgw)
H+	1.95618 × 10−6
Ca2+	0.30308
Mg2+	0.04185
Na+	0.3234
Cl−	0.551
HCO3−	0.007482

.

• Reactions with primary ((5)–(10)) and secondary minerals ((11) and (12)), using the Transition State Theory (TST)-derived rate laws [14]:

Calcite [CaCO_{3 (s)}] + H⁺_(aq) = Ca²⁺_(aq) + HCO₃⁻_(aq)

(5)

Quartz [SiO_{2 (s)}] = SiO_{2 (aq)}

(6)

Kaolinite [Al₂Si₂O₅(OH)_{4 (s)}] + 6.0 H⁺_(aq) = 5.0 H₂O_(aq) + 2 Al³⁺_(aq) + 2 SiO_{2 (aq)}

(7)

Albite [NaAlSi₃O_{8 (s)}] + 4 H⁺_(aq) = 3 SiO_{2 (aq)} + Al³⁺_(aq) + Na⁺_(aq) + 2 H₂O_(aq)

(8)

K-feldspar [KAlSi₃O_{8 (s)}] + 4 H⁺_(aq) = 3 SiO_{2 (aq)} + Al³⁺_(aq) + K⁺_(aq) + 2 H₂O_(aq)

(9)

Illite [K_0.6Mg_0.25Al_2.3Si_3.5O₁₀(OH)_{2 (s)}] + 11.2 H₂O_(aq) = 3.5 H₄SiO_4(aq) + 2.3 Al(OH)₄⁻_(aq) + 0.6 K⁺_(aq) + 0.25 Mg²⁺_(aq) + 1.2 H⁺_(aq)

(10)

Dolomite [CaMg(CO₃)_{2 (s)}] = 2 CO₃²⁻_(aq) + Mg²⁺_(aq) + Ca²⁺_(aq)

(11)

Siderite [FeCO_{3 (s)}] + H⁺_(aq) = HCO₃⁻_(aq) + Fe⁺⁺_(aq)

(12)

The mineral reactions are regulated by four parameters, as follows:

• o K_eq is the chemical equilibrium constant;
• o A is the reactive surface area;
• o E_a is the activation energy;
• o k₂₅ is the rate constant at T₂₅ = 298.15 K (25 °C), and the rate constant (k) in different temperature conditions can be expressed as follows:

$$k=k_{25}\times exp\left[-\frac{E_a}{R}\left(\frac{1}{T}-\frac{1}{T_{25}}\right)\right]$$

(13)

where R is the universal gas constant.

These parameters compose the TST dissolution and precipitation reaction rate (r), as follows [43]:

$$r=sgn\left[1-\left(\frac{Q}{K_{eq}}\right)\right]A{S}_w\left[k+{\sum }_{i=1}^{nct}{k}_i{a}_i^{wi}\right]{\left|1-{\left(\frac{Q}{K_{eq}}\right)}^{\xi }\right|}^{\zeta }$$

(14)

where

sgn: mathematical operator to return the sign of the expression;

Q: ion activity product (IAP);

S_w: aqueous phase saturation;

nct: number of reactant components;

a_i: activity of component i computed with Pitzer’s model [44,45];

wi: activity power;

$\xi ,\zeta$: TST model parameters.

The ratio log(Q/K_eq) is called the saturation index of the reaction. In this study, the numerical geochemistry model used assumes that if log(Q/K_eq) > 1, mineral precipitation occurs, and if log(Q/K_eq) < 1, mineral dissolution occurs. In the model used, the convention adopted is as follows: rate r is negative for dissolution and positive for precipitation. The rate is zero when log(Q/K_eq) = 1. Q and K_eq were assumed to be known from models [44,45] and open-source databases [46,47], respectively. However, A, k₂₅, and E_a were further investigated. The published values of the reactive surface area, for example, vary considerably among different authors due to the dependence on the grain size diameter of the rock minerals [17,48].

Table 3. Kinetic parameters of TST model for mineral reactions.

Mineral	Parameter	Value/Range	Source
Calcite	log10(Keq)	−8.38	[5,8,46,47,48,49]
	A	88 to 23,000 m2/m3
	log10(k25)	−8.8 to −0.3 mol/m2s
	Ea	14,400 to 41,870 J/mol
Quartz	log10(Keq)	−4.38	[5,8,46,47,49]
	A	2650 to 7128 m2/m3
	log10(k25)	−15 to −10 mol/m2s
	Ea	68,175 to 113,625 J/mol
Kaolinite	log10(Keq)	9.80	[5,8,46,47,50]
	A	1760 m2/m3
	log10(k25)	−13 mol/m2s
	Ea	62,760 J/mol
Albite	log10(Keq)	−19.7	[5,8,16,46,47]
	A	2563.38 m2/m3
	log10(k25)	−10.16 mol/m2s
	Ea	65,000 J/mol
K-feldspar	log10(Keq)	−22.66	[5,8,46,47,49]
	A	176 m2/m3
	log10(k25)	−12 mol/m2s
	Ea	67,830 J/mol
Illite	log10(Keq)	−43.28	[5,8,16,46,47]
	A	41,888 m2/m3
	log10(k25)	−10.98 mol/m2s
	Ea	23,600 J/mol
Dolomite	log10(Keq)	−16.46	[5,8,46,47,48,49]
	A	88 to 23,000 m2/m3
	log10(k25)	−9.22 to −3.19 mol/m2s
	Ea	36,100 to 62,760 J/mol
Siderite	log10(Keq)	−10.72	[5,8,46,47]
	A	4046.67 m2/m3
	log10(k25)	−3.19 mol/m2s
	Ea	36,100 J/mol

provides a concise overview of the kinetic values derived from an extensive literature survey. This uncertainty analysis will be specifically focused on the most representative and highly reactive minerals, namely quartz, calcite, and dolomite. For these minerals, ranges of values considered in the stochastic simulations are included in Table 3.

• Permeability alteration due to mineral precipitation or dissolution was computed by applying the Kozeny–Carman equation with an exponent value of 3, as Zeidouni et al. [51] recommended:

$$k_k=k_n/rf$$

(15)

where

$$rf={rf}_n{\left(\frac{{\phi }_n}{{\phi }_k}\right)}^3{\left(\frac{{1-\phi }_k}{{1-\phi }_n}\right)}^2$$

(16)

where the resistance factor rf is modeled by the Kozeny–Carman equation or the power–law relationship, k_n and k_k refer to permeability at previous (n) and current (k) timesteps, respectively. The porosity, φ, in (16) is calculated as follows:

$$\phi =\left[1+c_f\left(p-p^{*}\right)\right]\left[{\phi }^{*}-{\sum }_{j=1}^n\left(\frac{N_j}{{\rho }_{m,j}}-\frac{N_j^0}{{\rho }_{m,j}}\right)\right]$$

(17)

where φ* is the reference porosity without mineral precipitation/dissolution, N_j is the total moles of mineral j per bulk volume at the current time, $N_j^0$is the total moles of mineral j per bulk volume at time 0, ρ_m,j is the mineral molar density, c_f is the rock compressibility, and p* is the reference pressure.

2.2. Geological Model

A synthetic 2D model of a saline aquifer was constructed, where the petrophysical properties were assumed to be homogeneous. This model was used to assess methodologies for integrating geochemical reactions and to specifically emphasize the impact of geochemistry, excluding the influence of geological heterogeneity. The model contains a vertical injector perforated in the last three bottom layers, injecting 0.5 metric tonne/day of CO₂ over five years (1% of pore volume injected). Figure 1 shows the model size, discretized in 100 × 100 × 20 gridblocks with sizes of 10 m × 10 m × 5 m. The model discretization has been previously evaluated [22,32,52] and found to be suitable for representing reactive transport.

Regarding volumetric fraction, the mineralogy in this study consists of the following main primary minerals [53], which are considered relevant for the reactions and are included in the geochemical model (

Table 4. Volumetric fraction of the main primary minerals in the rock matrix.

Mineral Composition	Vol. %
Quartz	1.0
Kaolinite	0.15
Calcite	68.68
K-feldspar	0.12
Albite	0.05
Porosity	30.0

).

The relative permeability and capillary pressure curves for CO₂ and brine in carbonates were obtained from Bennion and Bachu [54]. Figure 2 shows the drainage curves used in the simulation of the aquifer.

The residual CO₂ trapping due to the relative permeability hysteresis with the saturation changes was modeled with the maximum gas trapped (S_gt) converted to the Land’s constant (C) [55] in the two-phase Carlson’s model [56], as recommended by Jarrell et al. [57], according to Equation (18):

$$S_{gt}=\frac{S_{gmax}}{1+C{S}_{gmax}}$$

(18)

where S_{g max} is the maximum gas saturation.

Burnside and Naylor [58] obtained S_gt distributions for CO₂ and brine in sandstone, shale, and carbonates based on more than 30 published coreflood data. The mean value of S_gt = 0.25 for carbonates was assigned in this study, which is equivalent to a C = 2.4.

Other properties of this reservoir model are summarized in

Table 5. Summary of reservoir properties.

Property	Value
Initial pressure	11,800 kPa @ 1200 m
Temperature	80 °C
Horizontal permeability	100 mD
Kv/Kh ratio	0.1
Porosity	30%
Rock compressibility	5.8 × 10−7 kPa−1

.

3. Results and Discussion

The objective functions (OF) that were evaluated are as follows:

• II: the CO₂ injectivity index after five years of injection, computed according to Equation (19):

$$II=\frac{q_{CO2}}{(p_{bh}-{\overline{p}}_{aq})},$$

(19)

where II is the CO₂ injectivity index, $q_{CO2}$ is the bottom-hole CO₂ volumetric injection rate; p_bh is the bottom-hole pressure, and ${\overline{p}}_{aq}$ is the average aquifer pressure.

• Mineralization_1000yr: Mineral dissolution or precipitation after 1000 years of CO₂ redistribution;
• Calcite_1000yr: Calcite concentration change (in gmoles) after 1000 years of CO₂ redistribution;
• Dolomite_1000yr: Dolomite concentration change (in gmoles) after 1000 years of CO₂ redistribution;
• Quartz_1000yr: Quartz concentration change (in gmoles) after 1000 years of CO₂ redistribution.

Table 6. Stochastics estimates for the objective functions using SR and MCS, associated errors, and R² (adjusted and predicted).

Objective Function	Estimate	SR	MCS	Error SR and MCS (%)	Adjusted R2 (%)	Predicted R2 (%)
		600 Runs	600 Runs	600 Runs	600 Runs	600 Runs
“II” (m3/d.kPa)	P10	0.0025560	0.0025560	0.0	99.7	86.1
	P50	0.0025550	0.0025550	0.0
	P90	0.0025510	0.0025510	0.0
“Mineralization_1000yr” (gmole)	P10	−2.286 × 106	−2.346 × 106	2.6	91.3	90.5
	P50	−2.377 × 106	−2.420 × 106	1.8
	P90	−2.795 × 106	−2.892 × 106	3.5
“Calcite_1000yr” (gmole)	P10	−2.286 × 106	−2.346 × 106	2.6	91.3	90.5
	P50	−2.377 × 106	−2.420 × 106	1.8
	P90	−2.796 × 106	−2.892 × 106	3.4
“Dolomite_1000yr” (gmole)	P10	57.70	55.2	4.3	82.0	80.0
	P50	53.45	53.0	0.8
	P90	51.71	50.4	2.5
“Quartz_1000yr” (gmole)	P10	It was not possible to calibrate a response surface because no changes were observed in relation to the initial content set for quartz throughout the simulation period.
	P50
	P90

summarizes the simulation results (SR) from 65,000 Monte Carlo simulations (MCS) using quadratic response surfaces (in Appendix A) as proxy models to represent the SR for each objective function (OF) in terms of the three estimates: P10, P50, and P90. In this context, the P90 represents the estimate with a 90% probability that the OF will have higher values. In other words, it signifies the lower end of the range of possible values for the OF, with only a 10% chance of the OF falling below this estimate. On the other hand, the P10 represents the upper end of the range of possible values for the OF, with only a 10% chance of the OF exceeding this estimate. The two last columns in the table display the respective R-squared (R²) coefficient, which measures the reduction in the variability of the response achieved by using the regressor variables in the model. The adjusted R² provides an estimate of the quality of fit of a proxy model to the simulation results (600 training experiments). On the other hand, the predicted R² is a measure of the model’s ability to accurately predict outcomes for new, unseen data (10 verification experiments), according to the representation shown in Figure 3. Overall, the results obtained from the SR and MCS analyses, along with the R² coefficients, indicate that the proxy models built using quadratic response surfaces are reliable and provide satisfactory estimates for the objectives of this study.

The errors between the SR and MCS results were found to be lower than 5%, with an R² higher than 80%. When considering only the P50 estimate, the prediction error decreased to less than 2%. These error levels are commonly observed in reservoir simulation applications [26,29,59]. Therefore, they were considered sufficiently accurate for this study, considering the numerous variables involved, and they will be considered for the following analysis.

3.1. Parametric Analysis with Proxy Models for Calcite Dissolution

Calcite dissolution is the primary mineral reaction observed in this case. Consequently, the “Mineralization_1000yr” and “Calcite_1000yr” outcomes from Table 6 are practically identical. In order to better assess the impact of the volumetric percentage of calcite (%vol) on mineral dissolution, an additional parameter was introduced into the experimental design to account for different proportions of calcite, ranging from 29% to 69%. The porosity was maintained at 30%, and as the amount of calcite decreased, the amount of quartz increased at the same rate. This was conducted to simulate a transition from a carbonate-rich to a siliciclastic-rich rock, where calcite essentially acts as the rock cementation.

Figure 4 provides a summary of the sensitivity analysis using a Tornado Plot, showcasing the most influential parameter on calcite dissolution. This includes the percentage of calcite in the original rock matrix, in addition to the kinetic parameters presented in Table 3. The rate constant (k₂₅) emerges as the most crucial parameter as it controls the reaction velocity and exhibits an explicit dependence on the energy of activation (E_a), as expressed in Equation (13) and illustrated in the interaction effect bar in Figure 4. The reactive area (A) follows as the second most influential parameter. The wide ranges of uncertainty associated with the parameters k₂₅ and A not only exert significant control over calcite dissolution but also interact with other parameters, thereby enhancing their overall impact. The broad range of A can be attributed to its dependence on the grain size diameter of the reactive minerals present in the rock, further contributing to the variability observed. In essence, the impact of an individual parameter on the overall objective function (OF) is known as the main effect. However, interaction effects take into account combinations of model parameters that cannot be described solely by the first-order effects. This finding aligns with the results drawn by Azuddin [20] in her history-matching study involving lab experiments on carbonate cores. In order to fit the data with numerical simulations, a tuning procedure was required for the rate constant and reactive area, as demonstrated in our results, which highlight these parameters as the most influential. This type of sensitivity analysis can support history-matching studies by identifying the parameters that need to be adjusted in numerical simulations.

In this particular analysis, the percentage of calcite exhibited a minimal effect. However, as will be discussed later, the limited amount of injected CO₂ is more influential in limiting the dissolution than the amount of calcite present. It explains why the difference between the P10 and P90 estimates is small for the calcite molar change.

The proxy model effectively replicates the temporal evolution of calcite dissolution. The comparison between the simulation results (in red) of a representative model for the P50 estimate and the predictions generated by the proxy model (in blue) is shown in Figure 5. In the context of reservoir simulation or modeling, representative models refer to a collection of models created to encompass the range of potential conditions or scenarios within the reservoir. These models are designed to capture various uncertainties and factors that may influence subsurface behavior, such as parameter variations or geological properties. In this particular case, the representative models correspond to the P10, P50, and P90 estimates, with the P50 estimate typically considered for project evaluation and approval.

This comparison demonstrates a notable agreement between the two sets of data, corroborating the accuracy and reliability of the proxy model. This alignment is further supported by the information presented in Table 6.

Table 7. Properties of the P50 representative model.

Property	Value	Unit
Rate Constant (k25)—Calcite (log10)	−8.44	mol/m2s
Reactive area (A)—Calcite	1202.26	m2/m3
Energy of activation (Ea)—Calcite	28,047	J/mol
Volumetric percentage (%vol)—Calcite	34	%

summarizes the properties of this representative model. In this particular case, a calcite percentage of 34% was assigned, taking into account the total dissolution after 1000 years, which amounts to −2.377 × 10⁶ gmoles. It is important to note that this dissolution represents only 0.0004% of the initial calcite present in the rock matrix. This limited dissolution primarily occurs in the vicinity of the wellbore (see Section 3.3).

3.2. Parametric Analysis with Proxy Models for Dolomite Precipitation

In contrast to the Tornado plot for calcite dissolution (Figure 4), the dolomite Tornado plot in Figure 6 exhibits a strong interaction among the parameters, as indicated by the orange bars. This can be attributed to the fact that dolomite, in this case, is a secondary mineral formed as a result of the dissolution of calcite induced by the injection of CO₂. The carbonate rock acts as a buffer, preventing a significant decrease in pH. This, coupled with the increase in HCO₃⁻ ions resulting from both the dissolved CO₂ in the brine solution and calcite dissolution, leads to dolomite precipitation. As calcite dissolves, it releases Ca²⁺ ions into the reservoir brine, triggering dolomite precipitation and consuming Mg²⁺ ions from the brine solution. Importantly, the injection of CO₂ induces the dolomitization process alongside calcite dissolution. This finding aligns with laboratory studies [9], which suggest that elevated temperatures exceeding 60 °C (in our case, 80 °C) initiate dolomite precipitation. It also corroborates other numerical studies [3], demonstrating the intrinsic relationship between calcite dissolution and dolomite precipitation. As a result of the strong dependencies among multiple parameters, the proxy model for dolomite demonstrates a lower quality (as indicated by the R² value in Table 6) compared to the calcite model.

It is important to emphasize that the proxy models developed using the proposed methodology were able to capture this significant relationship between calcite dissolution and dolomite. Additionally, the models demonstrated that quartz had no impact on this dynamic due to its low reactivity, as well as the other minerals present in rock matrix with small fractions. By incorporating the relevant parameters and their interactions, our models successfully captured the complex dependencies between these two processes. This highlights the effectiveness and reliability of our approach in accurately representing the intricate interplay between calcite dissolution and dolomite precipitation.

According to Figure 6, the percentage of calcite in the rock has a positive relationship as a main effect on dolomite generation. This means that as the amount of calcite increases, more Ca²⁺ ions can be released into the solution, facilitating the precipitation of dolomite. This process is expected to be more efficient in alkaline pH conditions, where the CO₃²⁻ ion is abundant and can support the growth of dolomite. To incorporate this aspect into the model, a new experiment was conducted, introducing the initial brine pH as a parameter with a range from 5 to 8 [60]. Figure 7 depicts the influence of pH as a parameter, displaying a strong interaction effect as it tightly stimulates dolomite generation in conjunction with the other parameters.

Unlike dolomite precipitation, the influence of pH on the calcite reaction is not as pronounced. This is because the main mechanism for calcite is dissolution, which is more likely to occur in the acidic region generated around the injector well (pH about 4), regardless of the initial pH. The acidic conditions favor the dissolution of calcite, leading to the release of calcium ions (Ca²⁺) into the solution. Therefore, while pH plays a significant role in dolomite precipitation, its effect on the calcite reaction is less significant. This observation is another crucial piece of evidence obtained from our study, which can be utilized to inform the design of future laboratory experiments. By identifying the key parameters that strongly influence the system, our results provide guidance on which parameters should be carefully measured and considered in experimental setups.

3.3. Thorough Proxy Analysis

The Tornado plots for calcite and dolomite, depicting the influence of various parameters after 1000 years, reveal that the rate constant emerges as the most significant kinetic parameter. This parameter holds direct control over the reaction rate, especially in the context of calcite dissolution. However, it is important to note that the rate constant also exhibits a wide range of uncertainty in open databases, as Dethlefsen et al. [12] highlighted in their analysis of batch simulation results.

Building upon the foundational research on uncertainty assessment conducted by Ates et al. [11], who suggested that geochemical reactions become more significant at higher CO₂ injection rates, a new experimental design incorporated the CO₂ injection rate as a variable to better understand its relevance in comparison to the rate constant. As a result, the final uncertainty analysis includes the same set of kinetic parameters for calcite, dolomite, and quartz, along with the percentage of calcite in the initial rock and the initial brine pH. Additionally, the injection rate varied over five years, ranging from 0.1 to 2.5 metric tonnes/day.

Figure 8 presents the updated Tornado plots for the objective functions (OFs) “Calcite_1000yr” and “Dolomite_1000yr”. The findings support the observations made by Ates et al. [11], as the CO₂ injection rate emerges as the most influential parameter. This parameter plays a crucial role in fueling the calcite dissolution reaction, which, in turn, is closely linked to the growth of dolomite, as evidenced by the right Tornado plot in Figure 8. Consequently, a positive relationship is expected and observed between the injection rate and both calcite dissolution (with a much greater magnitude) and dolomite precipitation. These results further underscore the significance of the CO₂ injection rate in driving mineral reactions and highlight its impact on the carbonate system.

The direct impact of calcite dissolution can be quantified by monitoring the evolution of injectivity, which serves as an important parameter in laboratory experiments and field-scale operations. By applying Equations (15)–(17), the rock porosity and permeability can be updated as a result of geochemical reactions, leading to an increase in the injectivity index (II). Considering the significant uncertainty associated with calcite dissolution, injectivity uncertainty is expected to also increase proportionally to the extent of calcite dissolution, particularly at higher injection rates, as illustrated in the left Tornado plot in Figure 8 and the crossplot in Figure 9.

Table 8. The ratio between P10 and P90 estimates for the CO₂ injectivity index after 5 years of injection.

Numerical Experiment	P10/P90 Ratio for II
Section 3.1 and Section 3.2: Fixed CO2 injection rate of 0.5 metric tonne/day	1.002
Section 3.3: CO2 injection rate as a parameter, varying from 0.1 to 2.5 metric tonnes/day	3.670

provides an assessment of the uncertainty level, measured by the ratio of II in the P90 estimate to that in the P10 estimate. By comparing the first experiment (Section 4), where the CO₂ injection rate was fixed at 0.5 tonne/day, with the second experiment, which incorporates a broad range of CO₂ injection rates as a parameter, the table highlights the uncertainty increase associated with injectivity in the second experiment, where greater amounts of CO₂ are involved. This information can be valuable in field and laboratory tests, as it aids in calibrating the balance between the dissolution and precipitation of carbonate minerals, with injectivity serving as a monitoring parameter, in addition to the effluent composition.

It is important to acknowledge that precipitation reactions of other secondary minerals, such as halite scale in the dry-out zone, can occur and partially counteract the changes in porosity and permeability resulting from calcite dissolution. Although this specific effect was not evaluated in the present study, it has been extensively studied and documented in work conducted by Machado et al. [22,61].

The significant variation observed in the permeability changes among the different scenarios depicted in Figure 10 reflects the uncertainty quantified in Table 8 when the CO₂ injection is a parameter. Furthermore, the primary changes observed in the model are concentrated at the top of the injection interval and in the shallowest layers. This observation raises concerns regarding the integrity of the caprock, as discussed in [62]. These changes can be attributed to the upward migration of CO₂ within the reservoir, driven by its buoyancy, resulting in a localized concentration of effects near the injection zone.

The uncertainty highlighted in this analysis emphasizes the range of potential outcomes and the importance of carefully considering the uncertainties associated with alterations in injectivity induced by calcite dissolution. Understanding and monitoring these uncertainties is crucial for ensuring the long-term effectiveness and integrity of CO₂ storage projects.

4. Conclusions

In conclusion, despite the use of a specific and simple geological model, the uncertainty analysis conducted in this study allows for the expansion and generalization of the main conclusions, providing valuable insights into carbonate reactions and informing future research and experimentation in this field. They are summarized as follows:

• The workflow employed in this study successfully mapped the uncertainty associated with the amount of CO₂ mineralized in carbonate rocks. This is essential for defining mitigation measures to address the risk of leakage in storage targets containing reactive minerals within their matrix;
• Proxy models constructed using quadratic response surfaces effectively represented the objective functions used to evaluate CO₂ injectivity and mineralization over extended periods of CO₂ redistribution. These models performed well even with a relatively small number of runs, making them valuable decision-making tools;
• The development of proxy models allowed for an expanded parametric analysis and facilitated the examination of the relationships among kinetic parameters and other key properties, such as the CO₂ injection rate, brine pH, and the initial percentage of calcite in the rock matrix. The findings highlight the significance of calcite dissolution as the primary mechanism of mineralization and the interplay between pH, available CO₂, and dolomite and calcite kinetics for the formation of dolomite as a secondary mineral;
• The parametric assessment revealed that the rate constant emerges as the most significant kinetic parameter. This parameter holds direct control over the reaction rate, especially in the context of calcite dissolution. Therefore, the uncertainty on this parameter underscores the need for further research and data collection to improve our understanding and its characterization, as highlighted in our results;
• Both pH and the initial concentration of calcite exert significant control on dolomite precipitation, while the CO₂ injection rate affects both calcite dissolution and dolomite generation. No changes in quartz content were observed in the conducted numerical experiments;
• The CO₂ injectivity index has demonstrated its capability to reflect the uncertainty associated with carbonate mineral reactions. Therefore, tracking the injectivity index from field or lab tests provides a practical means of monitoring and quantifying the levels of uncertainty involved.

These conclusions provide valuable support for the design of laboratory experiments, specifically regarding the measurement of geochemical reactions in corefloods. Parameters such as calcite amount, initial brine pH, and injection rate can be carefully controlled and varied in these experiments to enhance our understanding of carbonate mineral reactions.

Furthermore, to improve the accuracy of existing databases for characterizing mineral reactions induced by CO₂ injection in carbonate rocks, it is highly beneficial to further explore key parameters such as the rate constant and reactive area of the involved minerals. This can be achieved through the history matching of laboratory experiments, where these parameters are refined and adjusted to better align with observed data. By iteratively refining these parameters during the history-matching process, the precision and reliability of the current databases can be significantly improved. This, in turn, enhances our ability to accurately characterize and predict mineral reactions in carbonate rocks induced by CO₂ injection.

The findings yield valuable insights into the feasibility of CCS in carbonate saline aquifers by modeling the geochemical reactions and mapping the associated risks. This information is crucial for developing and implementing effective CCS strategies in geological formations that are not typically prioritized for carbon storage, such as carbonate rocks.