The Interplay between Experimental Data and Uncertainty Analysis in Quantifying CO2 Trapping during Geological Carbon Storage

Article Open Access

The Interplay between Experimental Data and Uncertainty Analysis in Quantifying CO2 Trapping during Geological Carbon Storage

Author Information
PETROBRAS, Rio de Janeiro, RJ 20231-030, Brazil
Hildebrand Department of Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin, TX 78712, USA
Authors to whom correspondence should be addressed.
Citations: 1
Clean Energy and Sustainability 2024, 2 (1), 10001;

Received: 05 December 2023 Accepted: 10 January 2024 Published: 16 January 2024

Creative Commons

© 2024 by the authors; licensee SCIEPublish, SCISCAN co. Ltd. This article is an open access article distributed under the CC BY license (

ABSTRACT: Numerical simulation is a widely used tool for studying CO2 storage in porous media. It enables the representation of trapping mechanisms and CO2 retention capacity. The complexity of the involved physicochemical phenomena necessitates multiphase flow, accurate fluid and rock property representation, and their interactions. These include CO2 solubility, diffusion, relative permeabilities, capillary pressure hysteresis, and mineralization, all crucial in CO2 trapping during carbon storage simulations. Experimental data is essential to ensure accurate quantification. However, due to the extensive data required, modeling under uncertainty is often needed to assess parameter impacts on CO2 trapping and its interaction with geological properties like porosity and permeability. This work proposes a framework combining laboratory data and stochastic parameter distribution to map uncertainty in CO2 retention over time. Published data representing solubility, residual trapping, and mineral trapping are used to calibrate prediction models. Geological property variations, like porosity and permeability, are coupled to quantify uncertainty. Results from a saline sandstone aquifer model demonstrate significant variation in CO2 trapping, ranging from 17% (P10 estimate) to 56% (P90), emphasizing the importance of considering uncertainty in CO2 storage projects. Quadratic response surfaces and Monte Carlo simulations accurately capture this uncertainty, resulting in calibrated models with an R-squared coefficient above 80%. In summary, this work provides a practical and comprehensive framework for studying CO2 retention in porous media, addressing uncertainty through stochastic parameter distributions, and highlighting its importance in CO2 storage projects. 
Keywords: CO2 trapping mechanisms; CCS; Uncertainty analysis; Proxy models; Saline aquifer

1. Introduction

Scenarios considered by international entities such as the International Agency (IEA) point out the complexity of the actions needed to achieve the goals established for reducing emissions of Greenhouse Gases, in the coming decades, simultaneously with the desired transformation in energy generation sources. The report “World Energy Outlook 2022” [1] emphasizes that multiple technologies and energy sources will play an essential role in providing energy resources for the planet in a more sustainable scenario. Among the various options for reducing CO2 emissions, Carbon Capture and Storage (CCS) presents itself as a technology with significant potential to reduce CO2 emissions in a probable scenario of continued use of fossil fuels in the coming decades. A CCS project involves CO2 capture from high-emission industries and injecting it into geological formations such as saline aquifers and depleted hydrocarbon reservoirs. One of the most critical barriers to long-term and large-volume CO2 storage in geological formations is the proof of safe and reliable storage. For that, CO2 can take advantage of different trapping mechanisms in porous media [2,3,4,5,6]; for example, free CO2 migration is controlled by the structural and stratigraphic trapping exerted by the caprock during the short-term encompassing the injection time, known as a primary mechanism [7]. Part of that mobile CO2 will be dissolved in the water (solubility trapping) as time passes. This CO2 solubility trapping is more pronounced in low-salinity brine, high-pressure, and low-temperature conditions [8,9]. The rising plume changes CO2 saturation leading to an increase in capillary trapped volume due to hysteresis in relative permeability and capillarity [10]. Some CO2 can be trapped as minerals because of brine pH reduction by CO2 injection of the brine and the mineralogy of the rock [11], on a timescale of more than 100,000 years [5]. Therefore, those additional mechanisms increase storage security when the buoyant CO2 is immobile in the pore space or no longer exists as a free phase (generally, as super-critical CO2). Modeling CO2 retention in porous media depends on the reservoir rock properties, such as porosity and permeability, as well as intrinsic parameters related to each trapping mechanism. These parameters can be measured in lab experiments using reservoir cores and reservoir brine samples. This data helps reduce the uncertainty in the amount of CO2 that will be trapped in a CCS project. However, it is important to note that a complete and single suite of lab data is not always available, or it may not fully capture the system’s heterogeneity, resulting in inherent modeling uncertainty [12]. Therefore, evaluating storage uncertainty involves assessing associated risks by considering modeling parameter uncertainties and system variability, enabling better understanding and management of CO2 storage risks. To address the evaluation of uncertainties in the context of CCS, various authors have focused on assessing uncertainties related to geological property distributions [13,14], such as porosity and permeability. Some studies have gone further by incorporating intrinsic parameters of CO2 trapping in their evaluations [15,16,17,18]. For instance, Jammoul et al. [17] conducted a history-matching study to reproduce experimental data while considering uncertainties. Similarly, Likanapaisal et al. [18] investigated the influence of rock and fluid properties, as well as grid size, on the shape of the CO2 plume. In their work, they specifically examined uncertainties related to the geometry of the CO2 plume, simplifying the modeling of trapping mechanisms by adopting a black-oil model. In this study, a framework for uncertainty analysis in CO2 retention in porous media is proposed. It integrates petrophysical properties and comprehensive trapping modeling, utilizing data collected from published laboratory studies. By doing so, it is possible to accurately quantify both trapped and free CO2 volumes, enabling a thorough analysis at a field scale with simulations spanning hundreds of years (long-term assessment). This is in contrast to previous studies [17], which focused on short-term redistribution over a few days in a smaller-scale lab model. Our study also emphasizes the importance of conducting risk assessments in geological carbon storage projects. It takes a comprehensive approach to characterize CO2 trapping mechanisms, including solubility, residual, and mineral trapping. This is achieved by leveraging lab data to establish distributions and relationships between parameters and geological properties. This unique approach sets our study apart from other studies conducted by different authors [15,16].

2. Modeling Trapping Mechanisms and Their Intrinsic Uncertainties

2.1. Solubility Trapping The solubility of CO2 in brine is Influenced by pressure, temperature, and salinity [8]. Several models have been developed to describe this relationship, such as those based on the Peng-Robinson Equation-of-state (EoS) and Henry’s law, as seen in [19] and [20]. Other models, like the one proposed by Duan and Sun [9], consider the EoS and the chemical potential of CO2 in both the liquid and vapor phases. This study will use the Li and Nghiem model [19], which is available in the CMG-GEM compositional simulator [21]. This model has been calibrated using published experimental data and it calculates Henry’s constant based on Equation (1), which takes into account pressure and temperature. Additionally, the impact of salt on gas solubility in the aqueous phase is accounted for by the salting-out coefficient [22].
```latex \text{ln}(H_i)=\text{ln}(H_i^∗)+\frac{\overline{v_i}}{RT}(p-p^∗) ```
where Hi: Henry’s constant at current pressure (p) and temperature (T); Hi: 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 (CO2 in this work). The parameters Hi and $$\overline{v_i}$$ can be adjusted to better fit the experimental solubility data and address any discrepancies in the model predictions. In this study, two levels of solubility were considered: 58,000 ppm at 50 °C [23] and 232,000 ppm at 80 °C [9]. As shown in Figure 1, the default values of Hi and $$\overline{v_i}$$ in Li and Nghiem model provide good agreement with the experimental results for the lower salinity case of 58,000 ppm. However, Figure 2 indicates that a tuning step is necessary to calibrate this model with the laboratory data when the brine salinity is higher than 230,000 ppm. Table 1 provides a summary of the changes in the regression parameters for each salinity case. These variations will be employed to quantify and prioritize the impact of the uncertainty in this parameter on the estimation of trapped CO2. This is particularly significant in the absence of laboratory data to use for solubility calculations. The synthetic water compositions are represented in Table 2.
Figure 1. CO2 solubility in brine with 58,000 ppm at 50 °C (in blue dots [23]) and prediction with Li and Nghiem model with default parameters (green curve) and after regression (red curve).
Figure 2. CO2 solubility in brine with 230,000 ppm at 80 °C (in blue dots [9]) and prediction with Li and Nghiem model with default parameters (green curve) and after regression (red curve).
Table 1. Solubility parameter changes for each salinity case.
Table 2. Ionic composition of formation brines.
In relation to solubility trapping, the dissolution of CO2 into brine can be accelerated by molecular diffusion through natural convection, as observed in the study by Rezk et al. [24]. This phenomenon significantly increases the overall storage rate in the aquifer, as fresh brine is brought to the top by convection currents. To assess the importance of molecular diffusion on storage results, a uniformly distributed effective diffusion coefficient of supercritical CO2 in water ranging from 3.5 × 10−5 to 4.0 × 10−5 cm²/s, based on the work of Ahmadi et al. [25], is used. This distribution encompasses the uncertainty arising from the dependence of the diffusion coefficient on temperature, pressure, and CO2 phase alteration (from gas to supercritical phase). However, it should be noted that it does not account for the potential variation of diffusion coefficient with brine salinity [26]. This is primarily due to the challenge of obtaining a comprehensive laboratory dataset that considers all these parameters. 2.2. Residual Trapping The residual trapping of CO2 resulting from relative permeability - capillary pressure - saturation hysteresis, was modeled using the two-phase Carlson’s model [27] with the maximum trapped gas (Sgt) converted to the Land’s constant (C) [28] in the two-phase Carlson’s model [27], as recommended by Jarrell et al. [29], according to Equation (2):
```latex S_{gt}=\frac{S_{g\;max}}{1+CS_{g\;max}} ```
where Sg max is the maximum gas saturation. To capture the variability of Sgt in sandstone, Burnside and Naylor [10] gathered data from over 30 published CO2-brine coreflood experiments. Using the published data, a distribution of Sgt (as depicted in Figure 3) is established and employed to assess its impact on the amount of residual CO2 trapped. The study by Burnside and Naylor [10] also presents a relationship between the endpoint relative permeability of CO2 (Max CO2 kr) and rock (absolute) permeability as shown in Figure 4. To account for the variability in rock permeability within geological uncertainties, the Max CO2 kr will be updated according to the permeability assigned to each gridblock. This ensures that the relative permeability curves can be adjusted accordingly for each gridblock in the uncertain geological scenarios considered in the study. The residual water saturation (Swr) is a critical parameter for determining the drainage relative permeability curves (kr) in two-phase CO2 and brine flow. Burnside and Naylor [10] measured a wide range of values for Swr. A continuous uniform distribution ranging from 0.15 to 0.5 was assumed for the subsequent steps in this study. Additionally, the hysteresis model proposed by Carlson [27] was considered in this study to generate the imbibition kr - Pc - Sw curves. In order to ensure consistency between the CO2-brine capillary pressure curves (Pc), kr, and the assigned Sgt, a script was developed. This script builds the Pc curves based on the power-law equation proposed by Brooks and Corey [30]. Equations 3 and 4 were used for the drainage and imbibition cycles, respectively. For the imbibition curve, the equation was adapted to respect the Sgt value as an endpoint kr (Figure 5). This adaptation ensures that the Pc curves are consistent with the Sgt value.
```latex P_{c,d}=(Max\;P_c -P_c^{S_{g \;crit}})(\frac{1-S_{g\;crit}-S_w}{1-S_{g\;crit}-S_{wr} })^{exp,d}+Min\;P_{c,d} ```
```latex P_{c,i}=(Max\;P_c)(\frac{1-S_{gt}-S_w}{1-S_{gt}-S_{wr} })^{exp,i} ```
where Sg crit is the critical gas saturation.
Figure 3. Prior Sgt distribution built with data collected by Burnside and Naylor [10]. The vertical green lines in the graph indicate the range of variation for Sgt.
Figure 4. Relationship between CO2 kr endpoint and the rock permeability based on the data collected by [10].
Figure 5. Schematic of CO2 and brine Relative Permeability (kr CO2 and krw) and Capillary Pressure (Pc) curves for drainage (subscript d) and imbibition (subscript i) processes.
In the course of numerical simulations, the Pc curves are calculated using J-functions [31]. These curves are then denormalized according to the rock permeability and porosity from the respective uncertain geological scenario. The resulting drainage kr and drainage/imbibition Pc curves are shown in Figure 6.
Figure 6. Drainage relative permeability curves for CO2 (A) and brine (B) and drainage (C) and imbibition (D) capillary pressure curves, in kPa, for CO2 and brine.
The curves in Figure 6 take into account the assumed distribution for Sgt (as depicted in Figure 3), the dependency of endpoint CO2 kr on rock permeability (as illustrated in Figure 4), and the Pc derived from the J-function considering various permeability and porosity values. Ranges of exponents were assigned based on laboratory data for the generation of relative permeability [32] and capillary pressure [33,34] curves. The table summarizing these ranges will be presented in the next section (Table 5), and the resulting curves will be used in the uncertainty analysis step of the study. In the current geo-model under investigation, the dry-out effect resulting from water vaporization, particularly around the well, has not been taken into consideration. However, the impact of this dry-out effect on injectivity has been studied by Machado et al. [35]. To model the dry-out effect, it is necessary to extrapolate the Pc curves to lower water saturation values (<Swr). Melrose [36] suggested a linear trend of Pc with the ln(Sw) when water is the wetting phase. This trend can guide the extrapolation process. 2.2. Mineral Trapping In this case, calcite cementation within the sandstone pores is assumed to be 1.1% [37]. The calcite reaction with CO2 can result in its dissolution and/or further precipitation. This process is governed by reactions occurring at the mineral/brine interface as follows: • The acidic reactions leading to bicarbonate and carbonate ions are controlled by kinetic parameters obtained from the PHREEQC database [38,39]:
OH + H+ = H2O
CO2 + H2O = H+ + HCO3
CO32− + H+ = HCO3
• Reactions with primary minerals, using the Transition-State-Theory (TST)-derived rate laws [40]:
Calcite [CaCO3] + H+ = Ca2+ + HCO3
Considering the findings of Machado et al. [35,41], it has been determined that only the reaction involving calcite will be included in the model. This decision is based on the fact that the other low-reactivity primary minerals, namely Quartz, K-feldspar, and Plagioclase, do not exhibit significant changes within the evaluated timeframe. The calcite reaction is regulated by four parameters, with three of them having similar values based on data obtained from open databases (PHREEQC: [38,39]; MINTEQ: [42]; WOLERY: [43]). However, the reactive surface area (A) changed considerably in value among different authors due to its dependence on the grain size diameter of the reactive minerals in the rock [44,45], especially in a clastic rock, such as sandstone, where the reactive surface area of an individual mineral grain is dependent on its grain size. • Keq is the chemical equilibrium constant for the calcite reaction; • A is the reactive surface area for calcite; • k25 is the rate constant of the calcite reaction at 25 °C (reference); • Ea is the activation energy. Table 3 presents a summary of kinetic values obtained from the literature survey that are used in the uncertainty evaluation.
After presenting the trapping modeling, it is worth stating the main assumptions in this paper. • No geomechanical modeling; • Pure CO2 stream is injected, e.g., without impurities and free water content; • No dry-out effect modeling due to water vaporization with CO2 injection; • No rock wettability changes; • No dependence of the diffusion on brine salinity; • Simplified geochemical modeling, only considering the calcite dissolution reaction; • The focus on the technical subsurface phenomena related to CO2 retention in porous media. It does not include the assessment of risks and uncertainties associated with surface processes [46], process efficiency [47,48], or economic aspects [49].
Table 3. Kinetic parameters of TST model for calcite reaction.

3. Uncertain Geological Scenarios and Results

The final evaluation of the CO2 storage potential in the sandstone aquifer was performed using a model based on actual seismic and well data from three wells in the Campos Basin, located in the state of Rio de Janeiro, Brazil. The aquifer covers both onshore and offshore areas. The sand grains in the aquifer are primarily composed of quartz and feldspars, with a notable presence of granitic lithoclasts, which are mechanically formed and deposited fragments of rock derived from older rocks. The average porosity of the sandstone is approximately 20%, and it is primarily generated through secondary processes such as matrix contraction, grain dissolution, or grain fracturing. The average permeability of the aquifer is around 2 Darcy. To evaluate the CO2 storage potential, 75 different petrophysical properties (some examples are given in Figure 7) are considered. Different realizations of permeability and porosity are generated to capture the variability and uncertainty in these properties (Personal communication with PETROBRAS, 2023) and final CO2 retention. These simulations provide insights into the potential capacity and effectiveness of CO2 storage in the sandstone aquifer, considering the range of geological and trapping properties observed in the region. This information is crucial for assessing the feasibility and long-term stability of CCS projects in the Campos Basin, contributing to the overall understanding of CO2 storage potential in this area. The reservoir model has a pore volume of 150 billion m³. However, for the specific analysis in this work, a sector model from the offshore region with a smaller pore volume of approximately 160 million m³ is considered. The gridblock size is 100 m × 100 m horizontally and 5 m thick. The grid refinement technique based on the workflow proposed by Machado et al. [41] is used to refine the grid near the CO2 injection well. The refinement is implemented to decrease the gridblock size from 100 m to 25 m near the injector. This adjustment aimed to improve the accuracy of representing the injection process and the trapping of CO2, whether as an aqueous component or a residual phase in porous media, or as a free supercritical fluid. The injection was evaluated for 20 years, from 2025 to 2045, with a constant injection rate of 1.0 million metric tonnes per year. A vertical injector was located in the center of the model in the bottom layers of the aquifer. It is important to note that in the selected area of the field, there is no evidence of an existing caprock (a layer of impermeable rock that can act as a barrier to prevent CO2 leakage to out of the storage zone, such as a shallow freshwater aquifer). As a result, the main risk associated with this project is the buoyancy of CO2 and the potential for leakage. Therefore, in this case, without the primary entrapment, it is crucial to address and quantify the uncertainty in the amount of CO2 trapped by the secondary mechanisms (mainly, by solubility and hysteresis) in the aquifer to assess the feasibility and safety of the project. Other properties of this actual reservoir model are summarized in Table 4. Table 5 provides a summary of the uncertainty range considered for the trapping properties.
Table 4. Summary of reservoir properties.
A numerical experiment was conducted with 400 simulations using Latin Hypercube design [50]. These runs combined 75 heterogeneous realizations of permeability and porosity, the ratio of Kv/Kh distribution, and the trapping parameters listed in Table 5. The numerical experiment took approximately 1.5 days to complete the simulations, utilizing 20 parallel jobs on a computer cluster with 40 processors (Intel® Xeon® Gold 6230 CPU @ 2.10GHz, 104 GB RAM). Additionally, another set of 1000 runs was carried out to validate the response surface approach. These simulations required an additional 2 days to complete. The four objective functions (OF) that were evaluated are as follows: • II: the CO2 injectivity index after 20 years of injection; • Trapped_shutoff: percentage of the total CO2 injected that is trapped due to solubility, residual, and mineralization after the injection stops; • Trapped_30yr: percentage of the total CO2 injected that is trapped after 30 years of redistribution; • Trapped_100yr: percentage of the total CO2 injected that is trapped after 100 years of redistribution.
Table 5. Uncertainty range of the trapping properties.
Figure 7. Cross-sectional view with four examples of x-permeability (A) and porosity (B) scenarios. The injector's location is indicated by the blue box.

4. Discussion

Table 6 summarizes the simulation results (SR) and from 65,000 Monte Carlo simulations (MCS) using quadratic response surfaces [51] as proxy models to represent the SR for each objective function (OF). The last column in the table displays the respective R-squared (R2) coefficient, which measures the reduction in the variability of the response achieved by using the regressor variables in the model. It is worth noting that the maximum amount of CO2 retention is reached after 30 years of redistribution. As a result, the uncertainty level calculated based on the ratio of P90 and P10 estimates, remains relatively constant at around 4 (in the 400-run case) and 3.5 (in the 1000-run case). In this context, the P90 represents the optimistic estimate for CO2 storage, while the P10 represents the pessimistic estimate. The errors between the SR and MCS results are found to be lower than 30%, with an R2 higher than 80%, when considering 400 runs. These error levels are commonly observed in reservoir simulation applications, as mentioned in the study by Cremon et al. [52]. Therefore, they were considered sufficiently accurate for this study, considering the numerous variables involved. Overall, the results obtained from the SR and MCS analyses, along with the R2 coefficients, indicate that the proxy models built using quadratic response surfaces are reliable and provide satisfactory estimates for the objectives of this study.
Table 6. Stochastics estimates for the Objective Functions using SR and MCS, associated errors, P90/P10, and R2.
The results indicate that the proxy models accurately predict the P50 estimate for the amount of CO2 trapped and the injectivity index. However, there are significant discrepancies between the P10 and P90 estimates predicted by MCS and those obtained from the cumulative probability distribution derived from the SR in Figure 8. These differences are particularly notable after 30 and 100 years of CO2 redistribution, with 400 simulation runs. This discrepancy leads to a reduction in the R2 from 95% to 83% for proxy models. To address this issue and better represent cases with lower retention (close to the P10 estimate) and higher retention (close to the P90), additional simulation runs are performed. However, it requires a larger number of simulations in the order of one thousand to achieve satisfactory agreement with the results from other authors [52,53]. Nevertheless, it is worth noting that despite the need for a larger number of simulations, the P50 estimate is well predicted by the proxy model even with the limited number of 400 runs. The P50 estimate is typically considered for project evaluation and approval.
Figure 8. Cumulative probability functions from SR (solid lines) and MCS (dashed lines) for 400 (in blue) and 1000 runs (in orange).
Table 7 presents the distribution of CO2 within the porous media for each uncertain scenario after 100 years. These distributions are calculated using MCS with a 1000-run case, and their respective R2 values are also provided. The table reveals that in most scenarios, the majority of the CO2 remains in a free state as a supercritical fluid. This finding highlights the potential risk of CO2 leakage to the surface in CCS projects, particularly in geological scenarios where the caprock cannot effectively act as a barrier. In these scenarios, when there is a larger amount of free CO2 present (P10), its plume has the potential to rise to shallower depths. This increases the risk associated with operations in areas that lack a caprock. Figure 9 illustrates the vertical distribution of CO2 as a free phase in models that represent the three different estimates: P10 (minimum plume depth: 1569 m), P50 (2114 m), and P90 (2320 m). This information obtained from the uncertainty analysis is highly relevant as it highlights the potential risk of the CO2 plume (under free CO2 form) undesirably coming into contact with shallow aquifers or even the sea floor (or surface in onshore targets). It underscores the importance of considering and managing these risks in CO2 storage projects. The P90 case represents the lower-risk scenario, as more than 50% of the injected CO2 is trapped, as shown in Table 7. Consequently, the free CO2 plume remains at deeper depths (Figure 9). It is important to note that in the timeframe studied here, there is no significant CO2 mineralized.
Figure 9. Depth of CO2 plume rise from equivalent simulation models for P10 (in red), P50 (in green), and P90 (in blue) estimates.
Table 7. CO2 inventory after 100 years of redistribution in porous media (OF: “Trapped_100yr”).
Figure 10 depicts the evolution of the percentage of CO2 trapped over a span of 120 years. The green symbols in Figure 10 represent simulation results corresponding to the P50 estimate, while the vertical lines represent the range between the P10 and P90 equivalent models. These equivalent models are simulation models that yield results similar to those obtained using the proxy model estimation. The plot highlights that the uncertainty during the injection period is relatively small since the contributions from secondary trapping mechanisms of solubility, capillarity, and mineralization are insignificant. However, as the injection ceases, the importance of these mechanisms gradually increases as CO2 redistributes within the porous media, which is consistent with Bachu [7]. This observation indicates that the uncertainty in CO2 trapping becomes more pronounced and influential over time after the injection period concludes. In the P50 estimate, it is observed that a plateau is reached from year 50 onwards, indicating the maximum amount of CO2 trapped through solubility and hysteresis. This suggests that 30 years of CO2 redistribution was sufficient to achieve the maximum entrapment for the given system. However, it should be noted that this plateau may vary if more optimistic (P90) or pessimistic (P10) trapping properties are assumed. MCS can also be employed to investigate the most influential parameters for CO2 retention. The effect estimated in a Tornado plot quantifies how a parameter (and some modification in its value) can impact the objective function (OF). By definition, the effect of an individual parameter on the OF is referred to as a main (or first-order) effect. However, interaction effects consider the combinations of model parameters that cannot be described by the first-order effects only [54]. In this particular problem, interactions are important and integral aspects of the models, such as those between relative permeability (kr) and capillary pressure (Pc) parameters such as residual saturations and relative permeability endpoint and exponents (Sgt, Swr, Sg crit, ng, and Max CO2 kr, as well as the petrophysical properties (K, ф, and Kv/Kh) derived from the geological models. Figure 11 presents the Top-7 most significant parameters and their associated interactions (represented by orange bars), which add up to 100% for the “Trapped_100yr” OF. If a linear response surface is employed as a proxy model, the R2 value reduces from 86% to 66% (in the case of 1000 runs), indicating the importance of capturing the interaction effects represented by the quadratic surface. This finding aligns with the conclusions of Khanal and Shahriar [16], who employed advanced machine-learning techniques to construct proxy models. Although their models had higher R2 values than the results of this work, they considered only seven input parameters and made simplified assumptions, i.e., no uncertainty in solubility, no capillarity effects, and the use of a single set of relative permeability curves.
Figure 10. Evolution of the % CO2 trapped over 20 years of injection and 100 years of redistribution. Green dots are the P50 estimate, the vertical black lines represent the variation between P10 and P90.
Figure 11. Tornado plot for the Objective Function “Trapped_100 yr” generated with MCS calibrated with 1000 runs.
In Figure 11, the solubility parameter Hi did not appear as a significant property for the objective function “Trapped_100yr”. This can be attributed to the narrow range (lower uncertainty) of the solubility parameter in the 58,000 ppm case, where the model prediction aligns well with the laboratory data, as shown in Figure 1. However, if a higher salinity case (230,000 ppm) is used, where the difference between the initial prediction and the data is more significant (as depicted in Figure 2), Hi will have a main effect of 3.7% on the “Trapped_100yr” objective function, indicating a significant uncertainty in solubility. These observations emphasize the importance of evaluating the uncertainty in the solubility model calculations. Figure 12 displays the Tornado Plot for “II” at the time of injection shutoff, specifically regarding CO2 injectivity. These results support the findings of Machado et al. [35], which concluded that residual CO2 saturation does not significantly impact injectivity. Essentially, the injectivity of CO2 is influenced by its effective permeability, which is determined by the product of the maximum CO2 relative permeability (Max CO2 kr) and the absolute permeability in both the horizontal (K) and vertical directions. The vertical permeability plays a critical role in enhancing CO2 buoyancy, which in turn facilitates its redistribution within porous media and increases its saturation. This amplification of CO2 saturation ultimately leads to an enhanced mobility of CO2 within the porous media, positively impacting its injectivity. The significant influence of geological properties on injectivity emphasizes the importance of the proposed methodology, which integrates different geological scenarios. Furthermore, Figure 12 indicates that the reactive calcite area, denoted as A, does not have a substantial effect because of small percentage of calcite of 1.1% present in the reservoir rock. The steps followed in the proposed workflow can be summarized by referring to the flowchart presented in Figure 13. The flowchart illustrates the integration of geological scenarios and trapping mechanisms, supported by experimental measurements, to map the fate of CO2 in porous media. This comprehensive approach enables a more precise assessment of the potential impacts of geological and trapping properties on CO2 injectivity and entrapment. Ultimately, it facilitates informed decision-making in CCS projects.
Figure 12. Tornado plot for the OF “II” generated with MCS calibrated with 1000 runs.
Figure 13. Flowchart with main steps for uncertainty analysis in CCS projects.

5. Conclusions

In conclusion, the uncertainty analysis conducted in this study, despite being based on a specific target with a geological model under different scenarios, allows for the expansion and generalization of the proposed workflow. This analysis provides valuable insights for project evaluation, particularly regarding the main subsurface risks associated with CCS operations. The main conclusions of this study can be summarized as follows: Regarding the experimental data: • Laboratory-measured data for parameters that control CO2 retention in porous media should be used in order to guide the input distribution for uncertainty assessment. This helps to more accurately map their impacts on the resulting trapped CO2 over time. Regarding the uncertainty assessment on the CO2 trapping: • The workflow employed in this study successfully mapped the uncertainty associated with the amount of CO2 trapped. This mapping is crucial for determining appropriate mitigation measures to address the risk of leakage in storage targets where the caprock may not serve as an effective barrier. It enables the evaluation of the feasibility of CCS operations in such targets and aids in making informed decisions. • The results obtained from a geological model of an actual saline sandstone aquifer demonstrated significant variations in CO2 trapping, ranging from 17% for the P10 estimate to 56% for the P90 estimate. In the P50 estimate, the CO2 entrapment reached a plateau of approximately 29% after 30 years following the cessation of injection. This highlights the importance of considering uncertainty during the risk assessment and approval phases of CO2 storage projects. • Proxy models constructed using quadratic response surfaces were able to accurately represent the objective functions used to evaluate CO2 injectivity and retention over different time periods in the P50 estimate. These models performed well even with a relatively small number of runs (400), making them valuable tools for decision-making. • Although the proxy models for the P10 and P90 estimates can benefit from further improvements, obtaining a larger number of additional simulations (more than 600 cases) can be challenging, especially for larger models. • Tornado plots revealed a strong interplay between reservoir properties (such as permeability - K -, and the ratio of vertical to horizontal permeability - Kv/Kh) and dynamic petrophysical parameters (such as Sgt, Swr, Sg crit, ng, and Max CO2 kr) in determining the amount of CO2 trapped as an objective function. This underscores the importance of an integrated workflow for uncertainty and risk evaluation. Regarding the uncertainty assessment on the CO2 injectivity: • The low reactivity of minerals in the rock matrix led to a negligible impact of mineral trapping, primarily attributed to calcite dissolution. However, this influence was not significant enough to affect CO2 injectivity; • By coupling the methodology with geological considerations, a more comprehensive understanding of the injectivity process can be achieved. The results highlight the greater influence of geological properties, such as the Kv/Kh ratio and absolute permeability, compared to trapping parameters. For future research, the authors recommend and intend to expand the geochemical modeling in the uncertainty evaluation to apply to targets where the matrix contains more reactive minerals. In this regard, the modeling of water vaporization and halite scale should be considered for a more representative and comprehensive study.


A: reactive surface area, m2/m3; exp: exponent for the capillary pressure correlation, dimensionless; Ea: activation energy, J/mol; Hi: Henry’s constant at current pressure (p) and temperature (T), dimensionless; $$H_i^∗$$: Henry’s constant at reference pressure (p) and temperature (T), dimensionless; J: Leverett J-function, dimensionless; K: rock permeability, mD; Keq: chemical equilibrium constant, dimensionless; k25: rate constant at 25°C, mol/m2s; krl: relative permeability for phase l, dimensionless; ng: CO2 kr curve exponent in Brooks-Corey model, dimensionless; nw: water kr curve exponent in Brooks-Corey model, dimensionless; Pc: CO2-brine capillary pressure, kPa; R: universal gas constant, 8.314 kPa·L/mol·K; Sgt: residual gas saturation due to hysteresis, dimensionless (vol/vol); Sg crit: critical gas saturation, dimensionless (vol/vol); Sl: current fluid saturation, dimensionless (vol/vol); Swr: residual water saturation, dimensionless (vol/vol); $$\overline{v_i}$$: partial molar volume at infinite dilution, L/mol; Greek symbols ф: rock porosity; Subscripts d: drainage; i: imbibition; l: fluid (w – water or CO2); h: horizontal; v: vertical. Acronyms EoS: equation-of-state; CCS: carbon capture and storage; II: injectivity index; MCS: Monte Carlo simulation; OF: objective function; SR: simulation results.


The authors express their gratitude to PETROBRAS for granting permission to use data and models in this study. They would also like to acknowledge Karim Salaheddine for his valuable contributions in collecting geochemical data and engaging in insightful discussions; and Dr. Chowdhury Mamun for thorough English revision.

Author Contributions

Conceptualization, M.V.B.M., M.D. and K.S.; methodology, M.V.B.M.; software, M.V.B.M.; validation, M.V.B.M., M.D. and K.S.; formal analysis, M.V.B.M., M.D. and K.S.; investigation, M.V.B.M.; resources, M.V.B.M. and M.D.; data curation, M.V.B.M. and M.D.; writing—original draft preparation, M.V.B.M.; writing—review and editing, M.D. and K.S.; visualization, M.V.B.M.; supervision, M.D. and K.S.; project administration, M.D. and K.S. All authors have read and agreed to the published version of the manuscript.

Ethics Statement

Not applicable.

Informed Consent Statement

Not applicable.


This research received no external funding.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Birol DF. World Energy Outlook 2022; IEA: Paris, France 2022.
Nghiem L, Shrivastava V, Kohse B, Hassam M, Yang C. Simulation of Trapping Processes for CO2 Storage in Saline Aquifers. In Proceedings of the Canadian International Petroleum Conference, Calgary, Alberta, June 2009.
Han WS, McPherson BJ, Lichtner PC, Wang FP. Evaluation of Trapping Mechanisms in Geologic CO2 Sequestration: Case Study of SACROC Northern Platform, a 35-Year CO2 Injection Site. Am. J. Sci. 2010, 310, 282–324. [Google Scholar]
Delshad M, Kong X, Tavakoli R, Hosseini SA, Wheeler MF. Modeling and Simulation of Carbon Sequestration at Cranfield Incorporating New Physical Models. Int. J. Greenh. Gas Control 2013, 18, 463–473. [Google Scholar]
Rackley SA. Carbon Capture and Storage, 2nd ed.; Butterworth-Heinemann: Cambridge, MA, USA, 2017.
Machado MVB, Delshad M, Sepehrnoori K. Potential Benefits of Horizontal Wells for CO2 Injection to Enhance Storage Security and Reduce Leakage Risks. Appl. Sci. 2023, 13, 12830. [Google Scholar]
Bachu S. CO2 Storage in Geological Media: Role, Means, Status and Barriers to Deployment. Progress Energy Combust. Sci. 2008, 34, 254–273. [Google Scholar]
Dodds WS, Stutzman LF, Sollami BJ. Carbon Dioxide Solubility in Water. Ind. Eng. Chem. Chem. Eng. Data Series 1956, 1, 92–95. [Google Scholar]
Duan Z, Sun R. An Improved Model Calculating CO2 Solubility in Pure Water and Aqueous NaCl Solutions from 273 to 533 K and from 0 to 2000 Bar. Chem. Geol. 2003, 193, 257–271. [Google Scholar]
Burnside NM, Naylor M. Review and Implications of Relative Permeability of CO2/Brine Systems and Residual Trapping of CO2. Int. J. Greenh. Gas Control 2014, 23, 1–11. [Google Scholar]
Xu T, Yue G, Wang F, Liu N. Using Natural CO2 Reservoir to Constrain Geochemical Models for CO2 Geological Sequestration. Appl. Geochem. 2014, 43, 22–34. [Google Scholar]
Machado MVB. Modelagem Numérica de Reservatórios de Petróleo: Prática Integrada de Simulação; Petrobras: Rio de Janeiro, Brazil, 2023.
Tadjer A, Bratvold RB. Managing Uncertainty in Geological CO2 Storage Using Bayesian Evidential Learning. Energies 2021, 14, 1557. [Google Scholar]
Mahjour SK, Faroughi SA. Selecting Representative Geological Realizations to Model Subsurface CO2 Storage under Uncertainty. Int. J. Greenh. Gas Control 2023, 127, 103920. [Google Scholar]
Raza Y. Uncertainty Analysis of Capacity Estimates and Leakage Potential for Geologic Storage of Carbon Dioxide in Saline Aquifers. Master’s Thesis, Massachusetts Institute of Technology: Cambridge, MA, USA, 2006.
Khanal A, Shahriar MF. Physics-Based Proxy Modeling of CO2 Sequestration in Deep Saline Aquifers. Energies 2022, 15, 4350. [Google Scholar]
Jammoul M, Delshad M, Wheeler MF. Numerical Modeling of CO2 Storage: Applications to the FluidFlower Experimental Setup. Transp. Porous. Med. 2023. doi:10.1007/s11242-023-01996-4.
Likanapaisal P, Lun L, Krishnamurthy P, Kohli K. Understanding Subsurface Uncertainty for Carbon Storage in Saline Aquifers: PVT, SCAL, and Grid-Size Sensitivity. In Proceedings of the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, October 2023.
Li Y-K, Nghiem LX. Phase Equilibria of Oil, Gas and Water/Brine Mixtures from a Cubic Equation of State and Henry’s Law. Can. J. Chem. Eng. 1986, 64, 486–496. [Google Scholar]
Harvey AH, Prausnitz JM. Thermodynamics of High-Pressure Aqueous Systems Containing Gases and Salts. AIChE J. 1989, 35, 635–644. [Google Scholar]
CMG. GEM Compositional & Unconventional Simulator (version 2022.10); Windows; CMG: Calgary, AB, Canada, 2022.
Bakker RJ. Package FLUIDS 1. Computer Programs for Analysis of Fluid Inclusion Data and for Modelling Bulk Fluid Properties. Chem. Geol. 2003, 194, 3–23. [Google Scholar]
Yan W, Huang S, Stenby EH. Measurement and Modeling of CO2 Solubility in NaCl Brine and CO2–Saturated NaCl Brine Density. Int. J. Greenh. Gas Control 2011, 5, 1460–1477. [Google Scholar]
Rezk MG, Foroozesh J, Abdulrahman A, Gholinezhad J. CO2 Diffusion and Dispersion in Porous Media: Review of Advances in Experimental Measurements and Mathematical Models. Energy Fuels 2022, 36, 133–155. [Google Scholar]
Ahmadi H, Erfani H, Jamialahmadi M, Soulgani BS, Dinarvand N, Sharafi MS. Corrigendum to “Experimental Study and Modelling on Diffusion Coefficient of CO2 in Water” Fluid Phase Equilibria 523 (2020) 112,584. Fluid Phase Equilibria 2021, 529, 112869. [Google Scholar]
Perera PN, Deng H, Schuck PJ, Gilbert B. Diffusivity of Carbon Dioxide in Aqueous Solutions under Geologic Carbon Sequestration Conditions. J. Phys. Chem. B 2018, 122, 4566–4572. [Google Scholar]
Carlson FM. Simulation of Relative Permeability Hysteresis to the Nonwetting Phase. In Proceedings of the SPE Annual Technical Conference and Exhibition, San Antonio, TX, USA, 4 October 1981.
Land CS. Calculation of Imbibition Relative Permeability for Two- and Three-Phase Flow from Rock Properties. Soc. Petrol. Eng. J. 1968, 8, 149–156. [Google Scholar]
Jarrell PM. Practical Aspects of CO₂ Flooding; Society of Petroleum Engineers: Richardson, TX, USA, 2002.
Brooks RH, Corey AT. Properties of Porous Media Affecting Fluid Flow. J. Irrig. and Drain. Div. 1966, 92, 61–88. [Google Scholar]
Leverett MC. Capillary Behavior in Porous Solids. Trans. AIME 1941, 142, 152–169. [Google Scholar]
Bennion DB, Bachu S. Drainage and Imbibition Relative Permeability Relationships for Supercritical CO2/Brine and H2S/Brine Systems in Intergranular Sandstone, Carbonate, Shale, and Anhydrite Rocks. SPE Reserv. Eval. Eng. 2008, 11, 487–496. [Google Scholar]
Krevor SCM, Pini R, Zuo L, Benson SM. Relative Permeability and Trapping of CO2 and Water in Sandstone Rocks at Reservoir Conditions. Water Resour. Res. 2012, 48, 2011WR010859. [Google Scholar]
Jung J, Hu JW. Impact of Pressure and Brine Salinity on Capillary Pressure-Water Saturation Relations in Geological CO2 Sequestration. Adv. Condens. Matter Phys. 2016, 2016, 1–11. [Google Scholar]
Machado MVB, Delshad M, Sepehrnoori K. Injectivity Assessment for CCS Field-Scale Projects with Considerations of Salt Deposition, Mineral Dissolution, Fines Migration, Hydrate Formation, and Non-Darcy Flow. Fuel 2023, 353, 129148. [Google Scholar]
Melrose JC. Valid Capillary Pressure Data at Low Wetting-Phase Saturations. SPE Reserv. Eng. 1990, 5, 95–99. [Google Scholar]
Xiao Y, Xu T, Pruess K. The Effects of Gas-Fluid-Rock Interactions on CO2 Injection and Storage: Insights from Reactive Transport Modeling. Energy Procedia 2009, 1, 1783–1790. [Google Scholar]
Parkhurst DL, Thorstenson DC, Plummer LN. PHREEQE: A Computer Program for Geochemical Calculations; U.S. Geological Survey: Denver, CO, USA, 1980.
Parkhurst DL, Appelo CAJ. Description of Input and Examples for PHREEQC Version 3: A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations.; U.S. Geological Survey Techniques and Methods, Book 6; U.S. Geological Survey: Denver, CO, USA, 2013; p. 497.
Lasaga AC. Chemical Kinetics of Water-Rock Interactions. J. Geophys. Res. 1984, 89, 4009–4025. [Google Scholar]
Machado MVB, Delshad M, Sepehrnoori K. A Practical and Innovative Workflow to Support the Numerical Simulation of CO2 Storage in Large Field-Scale Models. SPE Reserv. Eval. Eng. 2023, 26, 1541–1552. [Google Scholar]
Allison JD, Brown DS, Novo-Gradac KJ. MINTEQA2/PRODEFA2: A Geochemical Assessment Model for Environmental Systems: Version 3.0 User’s Manual; EPA: Washington, DC, USA, 1991.
Wolery TJ. EQ3/6, a Software Package for Geochemical Modeling of Aqueous Systems: Package Overview and Installation Guide (Version 7.0); USDOE: Washington, DC, USA, 1992.
Luo S, Xu R, Jiang P. Effect of Reactive Surface Area of Minerals on Mineralization Trapping of CO2 in Saline Aquifers. Pet. Sci. 2012, 9, 400–407. [Google Scholar]
Bolourinejad P, Shoeibi Omrani P, Herber R. Effect of Reactive Surface Area of Minerals on Mineralization and Carbon Dioxide Trapping in a Depleted Gas Reservoir. Int. J. Greenh. Gas Control 2014, 21, 11–22. [Google Scholar]
Chen H, Wang C, Ye M. An Uncertainty Analysis of Subsidy for Carbon Capture and Storage (CCS) Retrofitting Investment in China’s Coal Power Plants Using a Real-Options Approach. J. Clean. Prod. 2016, 137, 200–212. [Google Scholar]
Koelbl BS, Van Den Broek MA, Faaij APC, Van Vuuren DP. Uncertainty in Carbon Capture and Storage (CCS) Deployment Projections: A Cross-Model Comparison Exercise. Clim. Change 2014, 123, 461–476. [Google Scholar]
Van Der Spek M, Fout T, Garcia M, Kuncheekanna VN, Matuszewski M, McCoy S, et al. Uncertainty Analysis in the Techno-Economic Assessment of CO2 Capture and Storage Technologies. Critical Review and Guidelines for Use. Int. J. Greenh. Gas Control 2020, 100, 103113. [Google Scholar]
Oda J, Akimoto K. An Analysis of CCS Investment under Uncertainty. Energy Procedia 2011, 4, 1997–2004. [Google Scholar]
Mckay MD, Beckman RJ, Conover WJ. A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code.  Technometrics 2000, 42, 55–61. [Google Scholar]
Ma YZ. Quantitative Geosciences: Data Analytics, Geostatistics, Reservoir Characterization and Modeling; Springer International Publishing: Cham, Switzerland, 2019.
Cremon MA, Christie MA, Gerritsen MG. Monte Carlo Simulation for Uncertainty Quantification in Reservoir Simulation: A Convergence Study. J. Petrol. Sci. Eng. 2020, 190, 107094. [Google Scholar]
Ballio F, Guadagnini A. Convergence Assessment of Numerical Monte Carlo Simulations in Groundwater Hydrology: TECHNICAL NOTE. Water Resour. Res. 2004, 40, W04603. [Google Scholar]
Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global Sensitivity Analysis. The Primer, 1st ed.; Wiley: Wiley: Hoboken, NJ, USA, 2007.