Influence of Reservoir Heterogeneity on Simultaneous Geothermal Energy Extraction and CO₂ Storage

Abstract

This numerical study delves into the dynamic interaction between reservoir heterogeneity and its impact on the dual objectives of geothermal energy extraction and CO₂ sequestration. Employing finite element models, this research scrutinizes the effects of variable porosity, permeability, and capillary entry pressures on fluid dynamics and thermal processes within geothermal systems. Key findings reveal that these heterogeneities significantly dictate fluid behavior and heat distribution, influencing the operational efficiency and environmental sustainability of geothermal–CO₂ storage operations. By integrating the nonlinear, temperature-dependent properties of fluids, simulations provide in-depth insights into the coupled fluid–thermal interactions that govern system performance. The outcomes offer a refined understanding of the complex interdependencies within heterogeneous reservoirs, underpinning the optimization of design and operational methodologies for co-optimized geothermal energy and CO₂ storage solutions. Ultimately, this research contributes to the advancement of sustainable energy technologies, highlighting further investigative pathways to bolster the efficiency and longevity of two-phase water–CO₂ geothermal systems.

1. Introduction

The quest for sustainable energy solutions and climate change mitigation has led to an increased focus on geothermal energy extraction and carbon capture and storage (CCS) technologies [1,2,3,4]. The efficiency of these processes is significantly influenced by the intricate interplay of fluid flow, heat transfer, and phase behavior in porous media. In this context, the present study aims to provide a comprehensive understanding of how reservoir heterogeneity affects the efficiency of geothermal energy extraction and CO₂ storage, with a focus on sedimentary reservoirs.

Heterogeneity in sedimentary reservoirs arises from variations in rock properties such as permeability, porosity, and capillary pressure. Heterogeneity plays a crucial role in dictating fluid flow patterns, heat transfer mechanisms, and phase behavior. In the literature, considerable attention has been devoted to the intricate phenomena of multiphase flow and solute transport within porous media [5,6,7]. In a more specialized context, a subset of researchers has focused on modeling the dynamics of the dissolution of CO₂ [8,9,10,11]. A relatively limited number of investigations have extended this analysis to consider the dynamics of free-phase CO₂ within evolving capillary transition zones (CTZ), as well as the transport of dissolved CO₂ in brine [12,13,14]. Quantitatively, Singh et al. [13,14] demonstrated that under varying thermal conditions, the onset time for dissolved CO₂ solutal fingering was less than 10 years when the reservoir permeability κ exceeded 10⁻¹³ m². Given these empirical and theoretical findings, it becomes imperative not to overlook dissolution-trapping mechanisms in the early stages of a simulation. Motivated by these observations, this study aims to explore the concomitant effects of local capillary trapping on CO₂ storage over extended time scales.

The influence of porous media heterogeneity on capillary trapping mechanisms and the phenomena of gravitational and viscous fingering have been extensively investigated, albeit primarily in the realms of oil and gas extraction [15,16,17] and immiscible contaminant transport in aquifers [18,19,20]. Kueper et al. [18] conducted a series of drainage experiments in heterogeneous porous media. It was discovered that non-wetting phase fluids were unable to penetrate specific regions characterized by elevated capillary entry pressures. This observation was further substantiated by Kueper and Frind [19]. The findings indicated a high sensitivity of non-wetting phase fluid flow to capillary heterogeneity. Dawe et al. [21] performed drainage and imbibition displacement experiments in lensed and layered porous media and suggested that the propagation of the fluid front was predominantly governed by capillary heterogeneity. Chaouche et al. [16] corroborated this observation, noting significant alterations in saturation profiles that were attributable to capillary heterogeneity. Dale et al. [22] revealed a dichotomy in the dominant forces governing fluid displacement: at low flow rates, capillary forces were predominant, while at high flow rates, viscous forces took precedence. Van Duijn et al. [17] further elucidated the role of heterogeneity by demonstrating that even after drainage by a wetting phase fluid, the non-wetting phase (specifically oil) remained entrapped in regions of high permeability. Collectively, these studies indicate the incorporation of capillary heterogeneity into predictive models for CO₂ trapping in storage reservoirs.

Multiple methodologies have been employed to elucidate the complex dynamics involved during CO₂ storage in a heterogeneous saline reservoir. Han et al. [23] utilized numerical simulations to investigate buoyancy-driven ScCO₂ flow in heterogeneous reservoirs. The study indicated that a higher variance in permeability results in a smaller residually trapped mass of CO₂. Saadatpoor et al. [24] showed that CO₂ prefers to migrate through high-permeability pathways in the absence of capillary pressure variability. When capillary heterogeneity is introduced, CO₂ flows through areas of low capillary entry pressure, creating local capillary traps. Suekane et al. [25] reported a higher residual gas saturation in smaller glass beads, attributing it to elevated capillary pressures. Krevor et al. [26] found that the capillary barriers in sandstone rock could trap 2–5 times more CO₂ compared to residual trapping during water imbibition. Deng et al. [27] examined the uncertainties in storage capacity, CO₂ injectivity, and leakage. Injectivity was found to be influenced by local reservoir permeability, while leakage was a function of the connectivity of high-permeability facies. Ren et al. [28] and Rabinovich et al. [29] focused on the rate-dependent nature of CO₂ trapping. Ren et al. reported that locally trapped CO₂ saturation was independent of the injection rate, while capillary channeling and buoyancy-dominated flow regimes were rate-dependent. Trevisan et al. [30,31] indicated that small-scale heterogeneity can trap approximately seven times more CO₂ compared to a homogeneous sample. Gershenzon et al. [32] and Dai et al. [33] identified key parameters governing capillary trapping, including reservoir permeability and capillary entry pressure. Zhang et al. [34] employed 3D modeling, while Zulqarnain et al. [35] and Al-Khdheeawi et al. [36] emphasized the role of local capillary trapping in limiting CO₂ leakage. The latter also considered wettability heterogeneity.

Singh et al. [37] investigated the impact of heterogeneity on deep CO₂ storage reservoirs and its impact on the overall trapping of supercritical CO₂ (ScCO₂) during injection and post-injection periods. The study by Ershadnia et al. [38] focused on the impact of subsurface heterogeneity on geological CO₂ sequestration. The study evaluated how bottom-hole pressure (BHP), CO₂ breakthrough times, and trapping efficiencies were influenced by various factors including non-isothermal injection, geochemical reactions, capillary pressure heterogeneity, geomechanical effects, and permeability enhancement near the injection site. In the article by Guo et al. [39], the focus was on the impact of pore-scale surface wettability heterogeneity on immiscible two-fluid displacement in porous media. It was reported that at the continuum scale, the capillary pressure–water saturation (P_c–S_w) curves proved to be largely similar to those under homogeneous wetting conditions. This suggests that while pore-scale heterogeneities are crucial, their aggregate impact may sometimes manifest in ways that are similar to homogeneous systems. Fang et al. [40] revealed that CO₂ dissolution efficiency is significantly influenced by the level of permeability discrepancy within stratified formations. Extensive Monte Carlo simulations by Hansen et al. [41], on the other hand, highlighted the coupled effects of permeability heterogeneity and background flow on CO₂ trapping. Both studies underscored the complexity of CO₂ plume behavior in subsurface environments, indicating that increased heterogeneity can lead to a lower average and greater variance in trapping rates, which presents a challenge for CCS design. These investigations collectively suggest that the accurate prediction of CO₂ sequestration rates requires careful consideration of subsurface heterogeneity and background flow.

An understanding of subsurface temperature distribution and the rock properties of sedimentary formations is fundamental for evaluating and optimizing geothermal energy extraction processes [42]. These systems have previously been classified into high-temperature volcano-geothermal systems and low-temperature sediment-hosted geothermal systems [43]. Fluid migration pathways are pivotal in enhancing the efficiency and sustainability of geothermal energy extraction. Kolawole et al. [44] elucidated the mechanical zonation of rock properties and the development of fluid migration pathways in sedimentary-hosted geothermal reservoirs. Fan et al. [45] highlighted the depth of the Curie isotherm and the resistivity of sedimentary caprock as significant factors affecting a geothermal system’s performance.

Reservoir heterogeneity has a substantial impact on the flow dynamics and thermal performance of geothermal systems that utilize CO₂ as a working fluid [46]. In an endeavor to harness geothermal energy from sedimentary reservoirs, a series of studies has progressively unveiled the critical factors influencing system performance, especially when employing CO₂ as a working fluid [47]. Initially, Buscheck et al. [48] posited a preference for horizontal wells over vertical ones for enhanced injectivity and productivity, albeit without comparative metrics to substantiate this claim. This study was closely followed by the work of Elliot et al. [49] and Buscheck et al. [50], who proposed wellbore optimization strategies such as the incorporation of baffles and hydraulic ridges to augment system performance. Adams et al. [51] established a foundational understanding of the thermosiphon phenomena intrinsic to sedimentary geothermal reservoirs. Following this, Cihan et al. [52], alongside Jayne et al. [53] and Wang et al. [54], accentuated the role of geological heterogeneity in these systems. They noted the facilitation of faster thermal front propagation in heterogeneous formations but underscored the absence of rigorous mathematical models linking geological heterogeneity to breakthrough times and heat production efficiency. Babaei [55] explored the effects of boundary conditions and wellbore patterns on system performance. In a more recent attempt to unravel the intricacies of CO₂-based geothermal systems, Pan et al. [56] examined the temperature-dependent thermosiphon behavior of CO₂. This was followed by a notable contribution from Chen et al. [57], elucidating how steep geological slopes could amplify CO₂ injectivity, the CO₂ fraction in produced fluids, and the heat extraction rate, albeit with a caveat of early cold CO₂ breakthrough. This study further reiterated the complex interplay between geological heterogeneity and system performance, marking a continued evolution in understanding while highlighting persisting complexities and gaps in the literature warranting further investigation. The study by Singh and Chaudhuri [58] highlighted the influence of temperature and pressure on the properties of CO₂, revealing that in lower-temperature reservoirs, the enhanced heat capacity of CO₂ could lead to improved energy extraction without compromising sequestration efficiency.

The existing literature on both CO₂ trapping mechanisms and sedimentary geothermal reservoirs underscores the complex interplay of geological heterogeneity and thermodynamic variables such as porosity, permeability, capillary pressure, and heat transfer dynamics. On the one hand, the literature on CO₂ trapping mechanisms emphasizes the critical role of reservoir heterogeneity, stating that a quantitative understanding of variables like porosity, permeability, and capillary pressure is essential for accurate predictions and the effective management of CO₂ storage and migration. On the other hand, the literature on sedimentary geothermal reservoirs, while providing valuable insights, is noted for its gaps in rigorous mathematical modeling, particularly regarding thermodynamics and operational parameters. These two strands of research share a common theme: the need for a more nuanced, quantitative understanding backed by rigorous mathematical models and available data. Both fields would benefit from an integrated approach that combines rigorous mathematical modeling with empirical studies to address the complex interdependencies among geological, thermodynamic, and operational variables. Such an approach would not only strengthen the foundational understanding required for scaling and optimizing each respective system but could also reveal synergies and inform cross-disciplinary solutions for sustainable energy production and carbon capture and storage.

2. Methodology

For this study, we considered a geothermal system that was initially fully saturated with water and then tried to inject CO₂ from a wellbore opening that is located at the bottom of the system (Figure 1). Considering the buoyant flow of the CO₂, the production wellbore opening is considered to be at the top of the system. The same injection and production rates are used to constrain the operational plan. The study employs a robust mathematical framework to capture the dynamics of fluid flow in a two-phase water–CO₂ system within a heterogeneous reservoir [59]. A pivotal equation in this framework is the mass conservation equation (Equation (1)), which accounts for transient changes, advective transport, and dispersive effects. This equation is coupled with Darcy’s law (Equation (2)) to describe the fluid flow. Darcy’s law incorporates the effects of pressure gradient, gravity, and fluid properties, thereby providing a comprehensive description of the fluid dynamics in a porous medium. The coupling of these equations allows for a nuanced understanding of how reservoir heterogeneity influences fluid flow patterns, which is critical for both geothermal energy extraction and CO₂ storage.

$$\frac{\partial {\epsilon }_p{S}_f{\rho }_f}{\partial t}+\nabla .\left(S_f{\rho }_{f}u\right)=\nabla .(D_c\nabla (S_f{\rho }_f))$$

(1)

$$u=-\frac{k{k}_{r,f}}{{\mu }_f}({\nabla P}_f-{\rho }_{f}g)$$

(2)

To accurately model the capillary effects and relative permeability in a two-phase system, the study adopts the Brooks–Corey equations. Equations (3)–(6) provide a robust framework for describing the phase behavior in porous media. The Brooks–Corey model is particularly well-suited for capturing the pore-scale physics in heterogeneous reservoirs, making it an ideal choice for this study. These equations are integrated into the overall model to examine how capillary pressure and relative permeability contribute to the efficiency metrics of geothermal energy production and CO₂ storage capacity.

$$P_{{CO}_2}=P_w+P_c$$

(3)

$${\overline{S}}_n=\frac{\left(S_n-S_{rn}\right)}{\left(1-S_{rn}-S_{rw}\right)}; {\overline{S}}_w=\frac{\left(S_w-S_{rw}\right)}{\left(1-S_{rn}-S_{rw}\right)}$$

(4)

$$P_{S_n}=P_{S_w}+P_c\left(S_w\right); P_c\left(S_w\right)=P_{ec}{\overline{S}}_w^{-\frac{1}{{\lambda }_p}}$$

(5)

$$k_{{rS}_w}={\overline{S}}_w^{\left(3+\frac{2}{{\lambda }_p}\right)}; k_{{rS}_n}={\overline{S}}_n^2\left(1-{\left(1-{\overline{S}}_n\right)}^{1+\frac{2}{{\lambda }_p}}\right)$$

(6)

The energy conservation equation (Equation (7)) is employed to capture the heat transfer mechanisms in the system. This equation is comprehensive, accounting for transient thermal changes, advective heat transport, and various heat sources or sinks. To accurately model the effective thermal properties of the porous medium, the study uses these equations. The equations consider contributions from the fluid phase, the solid phase, and the immobile fluid phase, thereby providing a holistic understanding of heat transfer in heterogeneous reservoirs.

$${\left(\rho C_p\right)}_{eff}\frac{\partial T}{\partial t}+{\rho }_f{C}_{p,f}u.\nabla T+\nabla .q=Q+Q_p+Q_{vd}+Q_{geo}$$

(7)

$$q=-k_{eff}\nabla T$$

(8)

$${\left(\rho C_p\right)}_{eff}={\epsilon }_p{\rho }_f{C}_{p,f}+{\theta }_s{\rho }_s{C}_{p,s}+{\theta }_{imf}{\rho }_{imf}{C}_{p,imf}$$

(9)

$$k_{eff}={\epsilon }_p{k}_f+{\theta }_s{k}_s+{\theta }_{imf}{k}_{imf}+k_{disp}$$

(10)

The thermophysical properties of water, namely, dynamic viscosity (μ), specific heat capacity (C_p), density ($\rho$), and thermal conductivity (κ), exhibit temperature-dependent behavior that can be described by specific governing equations.

For dynamic viscosity, two distinct equations are applicable over different temperature ranges. Between 273.15 K and 413.15 K, the equation is:

$$\begin{gathered} \mu =1.38-2.12\times {10}^{-2}\times T^1+1.36\times {10}^{-4}\times T^2-4.65\times {10}^{-7}\times T^3+ \\ 8.90\times {10}^{-10}\times T^4-9.08\times {10}^{-13}\times T^5+3.85\times {10}^{-16}\times T^6 \left(273.15-413.15\mathrm{K}\right) \end{gathered}$$

(11)

For temperatures between 413.15 K and 553.15 K, the equation is:

$$\mu =4.01\times {10}^{-3}-2.11\times {10}^{-5}\times T^1+3.86\times {10}^{-8}\times T^2-2.40\times {10}^{-11}\times T^3 \left(413.15-553.15\mathrm{K}\right)$$

(12)

The specific heat capacity is governed by:

$${ C}_p=1.20\times {10}^4-8.04\times {10}^1\times T^1+3.10\times {10}^{-1}\times T^2-5.38\times {10}^{-4}\times T^3+3.63\times {10}^{-7}\times T^4$$

(13)

Density is described by:

$$\rho =1.03\times {10}^{-5}\times T^3-1.34\times {10}^{-2}\times T^2+4.97\times T+4.32\times {10}^2$$

(14)

Finally, thermal conductivity is governed by:

$$\kappa =-8.69\times {10}^{-1}+8.95\times {10}^{-3}\times T^1-1.58\times {10}^{-5}\times T^2+7.98\times {10}^{-9}\times T^3$$

(15)

Figure 2, sourced from the NIST website [60], provides a comprehensive depiction of the thermophysical properties of carbon dioxide (CO₂) as functions of both temperature and pressure. Specifically, Figure 2a illustrates the density ($\rho$), Figure 2b presents the dynamic viscosity (μ), Figure 2c shows the specific heat capacity (C_p), and Figure 2d displays the thermal conductivity (κ). These graphs offer valuable insights into how these properties vary under different thermal and pressure conditions. It is critical to note that unlike materials with constant or near-constant thermophysical properties, CO₂ demonstrates substantial variations in $\rho$, μ, C_p, and κ over a range of temperatures and pressures.

In the context of our study, a set of key parameters have been defined to construct the reservoir model, and these are summarized in

Table 1. Parameters required for modeling.

Parameter	Value	Parameter	Value
Injection rate	Case dependent	Mean porosity	0.2
Injection temperature	Case dependent	Mean permeability	5 mD
Mean capillary entry pressure	30 kPa	Residual CO2 saturation	0
Residual water saturation	0.2	Rock density	2600 kg/m3
Thermophysical properties of CO2 and water	f(p, T), NIST webbook	Thermal gradient	30 K/km
Rock thermal conductivity	3 W/m/K	Rock-specific heat capacity	850 J/kg/K
Wellbore radius	0.2 m	Wellbore length	50 m
Surface temperature	10 °C	System thickness	300 m
System length	500 m	System width	500 m

. These parameters serve as the foundational elements for the numerical simulations, ensuring a rigorous and quantitatively defined modeling approach. It should be noted that these data pertain to an imaginary reservoir, which serves as a representative model for real-world scenarios.

In our study, geostatistical modeling was performed using Stanford Geostatistical Modeling Software (SGeMS), a software suite primarily designed for geostatistical applications [61]. This software was chosen for its robustness and flexibility in handling various geostatistical simulation algorithms. We employed a sequential Gaussian simulation (SGS) algorithm. SGS is widely used for its efficiency in modeling continuous variables and its ability to reproduce spatial variability [62]. The algorithm involves transforming the data into a standard Gaussian distribution, performing simulations at unsampled locations, and back-transforming the results to the original data distribution. The variogram model used in this study was the spherical model. The spherical model is characterized by a gradual increase in spatial variance with distance until a defined range is reached, beyond which the variance stabilizes [63]. This model is chosen for its suitability for representing the spatial correlation of the studied geological attributes, providing a balance between short-range and long-range spatial dependencies. The model resolution, defined by the cell size in the simulated grid, was set to 25 m, 35 m, and 45 m. This resolution was determined based on a compromise between computational efficiency and the level of detail required for the accurate representation of spatial variability. Furthermore, a statistical seed of 1 was used for the random number generator in SGeMS, ensuring the reproducibility of the simulation results. SGeMS generates random numbers with an average of 0 and a variance of 1. These numbers are used through $k_i=k_{ave}\mathrm{e}\mathrm{x}\mathrm{p}(\delta \times n_i)$ to the assigned permeability values for the simulation nodes, where $k_{ave}$, $\delta$, and $n_i$ show the average permeability of the system, interested variance for the permeability field, and resulting number from SGeMS, respectively.

A hydrothermal model for reservoir-scale analysis was constructed using COMSOL Multiphysics software, version 6.1 [59]. The domain was discretized using free tetrahedral elements, and the mesh size was carefully chosen to strike a balance between result accuracy and computational speed. Accordingly, the maximum element size, minimum element size, maximum element growth rate, curvature factor, and resolution of narrow regions are set to 20 m, 9 m, 1.5, 0.6, and 0.5, respectively. With these specifications, the mesh effectively captures the intricacies of the model while maintaining computational efficiency. The model domain consists of 159,094 elements, while the boundary and edges contain 6612 and 266 elements, respectively. For the sake of numerical accuracy, the model adopts a relative tolerance of 0.001 and a tolerance factor of 0.1. The backward differentiation formula (BDF) has been selected for numerical discretization, and automatic time stepping is integrated into the computational process. An implicit solving approach has been employed. It should be noted that this hydrothermal framework has previously been validated for geothermal applications with water as the working fluid, as reported by Mahmoodpour et al. [64]. Moreover, the effectiveness of the implemented methodology and tool in simulating the CO₂-water system has been previously validated in the documentation of COMSOL Multiphysics^® under the title “CO₂ Storage in a Geologic Formation.” This documentation delves into the examination of CO₂ propagation in the Johansen formation, off the coast of Norway. The methodological setup provides a robust, quantitatively defined avenue for simulating hydrothermal processes at the reservoir scale. The numerical parameters, such as element sizes and tolerance values, have been meticulously selected to balance computational speed with the necessary resolution for capturing essential hydrological and thermal attributes.

Considering the data presented in Table 1, the base case is characterized by an injection and production flow rate of 15 kg/s, along with an injection temperature of 320 K. To assess the influence of flow rate and permeability, two cases have been established by adjusting the flow rate to 30 kg/s and the permeability to 2.5 mD. To explore the impact of heterogeneity, two variance levels (0.25 and 0.5) are examined across three correlation lengths: 25 m, 35 m, and 45 m. Simulations for each combination of variance and correlation length are repeated for 5 realizations.

3. Results

In Figure 3, we present the stochastic heterogeneous fields for permeability nodes, as generated using Stanford Geostatistical Modeling Software (SGeMS), with a focus on varying correlation lengths and a fixed variance of 0.25. The figure is organized into a matrix format, featuring three columns and five rows. Each column represents one of three distinct correlation lengths—25 m, 35 m, and 45 m. Conversely, the five rows illustrate five different realizations for each of these specified correlation lengths. Five realizations are performed to understand the averaged behavior of the heterogeneity effect. A careful examination of Figure 3 reveals distinct spatial patterns of permeability that are highly dependent on the correlation length. As the correlation length increases from 25 m to 45 m, there is a noticeable trend toward larger connected regions of similar permeability values, thereby implying a higher degree of spatial continuity. This is particularly evident when comparing the first column (correlation length of 25 m) with the third column (correlation length of 45 m). At the shorter correlation length of 25 m, the permeability fields exhibit higher local variability and less continuity, which manifests as a patchier distribution. Conversely, at a correlation length of 45 m, the permeability field shows more extended zones of similar permeability, highlighting the role of correlation length in influencing spatial heterogeneity. These realizations provide a quantifiable measure of the uncertainty inherent in the permeability distribution, which can be used to perform two-phase water-CO₂ fluid flow modeling during geothermal energy extraction. Similar to Figure 3, Figure 4 is structured into three columns and five rows. Each column corresponds to a unique correlation length—25 m, 35 m, and 45 m—while each row portrays a different realization at these specific correlation lengths. Upon scrutinizing Figure 4, it is evident that an increase in variance from 0.25 (as seen in Figure 3) to 0.50 leads to a broader spread of permeability values within each realization. This manifests as more pronounced contrasts in the spatial distribution of permeability, thereby implicating a heightened level of uncertainty in operations and difficulties in reservoir modeling, be it for CO₂ trapping mechanisms or sedimentary geothermal reservoirs. The influence of increased variance is consistent across all correlation lengths, but its impact is particularly noticeable at shorter correlation lengths where the system already exhibits high local variability. The relationship between variance and correlation length in Figure 4 resonates with the key objective of this study by illustrating how geological and thermodynamic variables are interdependent. For example, the elevated variance could potentially affect capillary trapping efficiency in CO₂ storage or the heat transfer dynamics in geothermal reservoirs. Similarly, the correlation length influences the spatial distribution of permeability, which, in turn, affects porosity and capillary pressure—key variables for CO₂ trapping and geothermal energy extraction. These processes will be discussed later in this work.

3.1. Influence of the Injection Rate

To elucidate the dynamics within a homogeneous reservoir, Figure 5 plays a critical role by portraying the distributions of supercritical CO₂ saturation, fluid temperature, and reservoir pressure over time for a base case with an injection rate of 15 kg/s and a uniform mean permeability of 5 mD. Figure 5, which is organized into three columns corresponding to the aforementioned variables and five rows representing time snapshots at 1, 2, 5, 15, and 30 years, offers a window into the thermodynamic processes at play in an idealized, uniform geological setting. As the simulation time progresses, the CO₂ plume is seen to spread outwards from the injection point, as shown by the left column, being influenced solely by the injection conditions and the intrinsic properties of CO₂ and the reservoir rock. The absence of heterogeneity simplifies the interpretation of CO₂ migration, which is primarily driven by pressure gradients and the phase behavior of CO₂ in response to the reservoir conditions. The middle column of the figure demonstrates how the temperature within the reservoir evolves over time in response to the introduction of supercritical CO₂. Initially, a disturbance in temperature fields is apparent, but as time advances, a quasi-steady state is approached, which aids in understanding the baseline heat transfer dynamics when devoid of the complexities introduced by geological irregularities. This insight is particularly crucial for the assessment of geothermal energy recovery in sedimentary basins, where homogeneous reservoir assumptions can serve as a first-order approximation for system behavior. The rightmost column, illustrating reservoir pressure, sheds light on the impact of CO₂ injection on the pressure regime of the reservoir. Over the 30-year period depicted, we can observe the trend toward pressure equilibration throughout the reservoir. This equilibration process, uninhibited by heterogeneous permeability barriers or preferential pathways, provides a fundamental perspective on the influence of injection operations on reservoir pressure—a factor that has implications for both the mechanical stability of the reservoir and the efficiency of CO₂ storage. Therefore, Figure 5 presents a baseline scenario of a homogeneous reservoir, offering valuable baseline data against which the complexities observed in heterogeneous systems can be compared. This comparison is vital for enhancing CO₂ sequestration strategies and for optimizing geothermal energy extraction in more complex geological formations.

In contrast to Figure 5, Figure 6 delves into a scenario with a doubled flow rate, where the injection rate is elevated to 30 kg/s while maintaining the same permeability. Figure 6, when juxtaposed with Figure 5, depicts a broader spread of supercritical CO₂ saturation across the reservoir. This enhanced dispersion, which can be attributed to the increased injection rate, suggests that the dynamics of CO₂ migration are markedly influenced by the rate at which CO₂ is introduced into the system. The elevated flow rate visibly amplifies the spread and uniformity of CO₂ saturation, underlining a direct relationship between the injection rate and the spatial extent of CO₂ plume development. As for fluid temperature, the findings from Figure 6 suggest a hastened approach toward thermal equilibrium within the reservoir, contrasting with the more gradual thermal diffusion observed in Figure 5. This rapid temperature stabilization, driven by the higher flow rate, emphasizes the dominant role of advection in heat transfer within the reservoir, a factor that gains significance in the design and management of sedimentary geothermal systems. Turning to reservoir pressure, the higher injection rate illustrated in Figure 6 is correlated with more pronounced pressure gradients, particularly when proximal to the injection site. This observation underscores the impact of injection rate on the mechanical behavior of the reservoir, necessitating refined pressure management strategies to uphold structural integrity and prevent deleterious geomechanical responses such as fracturing or caprock compromise. Temporally, the results from Figure 6 convey an expedited response of the reservoir system to the injection process. The quickened pace of this change, relative to Figure 5, necessitates a forward-looking approach to reservoir management that anticipates and adjusts for the more immediate impacts of higher injection rates on reservoir performance and longevity. This accelerated dynamic is a critical consideration for the strategic planning of CO₂ sequestration initiatives and the operational efficiency of geothermal energy systems.

3.2. Influence of Reservoir Permeability

Figure 7 scrutinizes the impact of halving the mean permeability to 2.5 mD while maintaining the same injection rate with respect to the base case (Figure 5). In this reduced permeability scenario, Figure 7 demonstrates a pronounced concentration of supercritical CO₂ within a more limited region of the reservoir. The lower permeability acts as a barrier to the lateral dispersion of CO₂, resulting in heightened localized saturation and potentially impacting the volumetric efficiency of CO₂ storage. This confinement suggests that the reservoir may require additional injection points to achieve the desired sequestration goals, highlighting the crucial role of permeability in designing storage strategies. Thermal behavior within the reservoir, as indicated by Figure 7, undergoes a more gradual transition to equilibrium. The limited permeability dampens the convective heat transfer, which, in turn, suggests a more protracted period of thermal propagation throughout the reservoir. This finding has ramifications for sedimentary geothermal reservoir operations, where the slower thermal response may necessitate a reevaluation of heat extraction methodologies and system design to ensure optimal performance. Regarding pressure dynamics, Figure 7 presents an increase in pressure gradients near the injection site when compared to Figure 5. The heightened resistance to fluid flow due to the reduced permeability necessitates a careful consideration of pressure buildup and its potential geomechanical consequences, such as the integrity of the caprock and the activation of faults. When observing the temporal progression of the aforementioned parameters in Figure 7, there is a discernible deceleration in the dynamic response of the reservoir. The permeability reduction imposes a lag on both CO₂ plume evolution and the thermal stabilization process, the implications of which extend to the operational planning of CO₂ sequestration and geothermal resource development. In essence, a comparative analysis of Figure 5 and Figure 7 provides insight into the influence of permeability on CO₂ distribution and thermal dynamics within a homogeneous reservoir system. Lower permeability delineates a scenario with constrained CO₂ migration, slower thermal equalization, escalated localized pressures, and overall decelerated system dynamics, all of which are pivotal considerations for the optimization and management of geothermal energy exploitation and carbon storage.

3.3. Local Capillary Trapping

Figure 8 presents an intricate portrayal of the local capillary trapping phenomenon that underscores the geological complexity and thermodynamic intricacies within the reservoir. The figure is organized into multiple panels that collectively offer a detailed analysis of the spatial and temporal dimensions of supercritical CO₂ trapping mechanisms. Figure 8a provides a detailed visualization of the permeability distribution within the first realization at a correlation length of 25 m and a variance of 0.50. This panel is pivotal as it lays the foundation for understanding the flow pathways and the spatial heterogeneity that influence CO₂ migration and entrapment. Figure 8b maps the capillary entry pressure across the reservoir, delineating the thresholds at which CO₂ invades the pore spaces. The heterogeneity in these pressures can lead to differential trapping across the reservoir, which is essential for evaluating the efficiency and safety of CO₂ sequestration projects. Subsequently, Figure 8(c1–f1,c2–f2) convey the evolution of CO₂ saturation and temperature distributions within the reservoir over time frames of 1, 2, 5, and 15 years. Red patches indicating stable zones of CO₂ saturation, attributable to local capillary trapping, are of particular interest. These patches, which remain consistent over time, exemplify areas where the CO₂ is sequestered by capillary forces, thereby mitigating its upward migration and contributing to the overall trapping capacity. The local capillary trapping process, as elucidated by these panels, has significant implications for the management of CO₂ storage projects. The stability of these trapped zones directly influences the long-term security and efficacy of CO₂ storage solutions [34,37]. Furthermore, this trapping mechanism bears significant implications for geothermal energy extraction. The presence of supercritical CO₂ pockets, shaped by local capillary forces, introduces heterogeneity to the thermal properties within the reservoir. Depending on the relative thermophysical properties of CO₂ and the geothermal brine, these trapped pockets can alter the thermal conductivity and heat transfer rates within the reservoir. The heterogeneity of the reservoir, as depicted by the permeability distribution (panel 8(a)) and the capillary entry pressure (panel 8(b)), directly influences the temperature evolution within the reservoir. This influence is due to the impact of heterogeneity on fluid flow and the subsequent distribution of heat. Areas with higher permeability facilitate greater fluid movement, which can lead to more efficient advective heat transfer, as the movement of fluids can transport heat away from the injection points more effectively. Conversely, areas with low permeability may act as thermal barriers, impeding fluid flow and, thus, heat transfer, leading to a more localized heating effect. Moreover, the heterogeneity of capillary entry pressure, which causes the local capillary trapping of supercritical CO₂, creates zones with varying saturation levels. These zones have different thermal conductivities and specific heat capacities, leading to an uneven thermal response when subjected to heat transfer processes [37]. The trapped CO₂ pockets may act as thermal insulators due to their lower thermal conductivity compared to the surrounding brine, potentially leading to a cooler region within the reservoir. This phenomenon necessitates a strategic re-evaluation of operational parameters to ensure the consistent and efficient extraction of geothermal energy. Adjustments may include modified pumping rates or the implementation of thermal stimulation techniques to counteract the reduced permeability and altered pressure gradients caused by CO₂ trapping. This dynamic interplay between CO₂ trapping and geothermal operation suggests the need for adaptable and informed reservoir management strategies to balance the objectives of CO₂ sequestration and energy extraction, potentially leading to increased operational complexities and costs.

In the quest for a more quantitatively rigorous understanding of geothermal reservoirs and CO₂ trapping mechanisms, Figure 9, Figure 10 and Figure 11 show the ratio of CO₂ (loosely termed as gas, here) to water content, calculated at the production well as a function of time. The gas-to-water ratio is a crucial parameter that can be linked to the effectiveness of heat extraction and CO₂ trapping. A lower ratio may indicate greater CO₂ storage but could also signify a reduced efficiency in geothermal energy extraction due to the reasons elucidated above. The variations in these plots for different correlation lengths and variances reflect the highly sensitive nature of geothermal and CO₂ storage systems regarding geological heterogeneity. This sensitivity validates the research objective, emphasizing the necessity for rigorous mathematical and empirical modeling to predict and manage the complex interdependencies in such systems. Furthermore, differences in the trajectories for the various realizations within each correlation length indicate that stochastic variability can have a substantial impact, identifying the need for probabilistic approaches in reservoir modeling and management.

Figure 12 provides a quantitative elucidation of the dynamic interplay between CO₂ and water in geothermal reservoirs over time, focusing on the gas-to-water ratio and the cumulative heat extracted at the production well. The changes in supercritical CO₂ saturation over time shown in Figure 5 can be directly related to the gas–water ratio curve in Figure 12. For instance, areas of high CO₂ saturation would correspond to higher values of GWR. The fluid temperature profiles in Figure 5 can be linked to the ‘cumulative heat extracted’ curve in Figure 12. High-temperature zones would contribute to a steeper slope in heat extraction lines. The pressure distribution in Figure 5 could influence the curves representing the contributions of water and CO₂ to the total extracted heat. High-pressure zones might facilitate more effective heat extraction from water, thereby affecting the heat extracted by individual fluids.

Figure 13 provides a multi-faceted, quantitative assessment of the interactions between CO₂ and water within a geothermal reservoir in the specific case of realization 1, with a correlation length of 25 m and a variance of 0.50. The figure encapsulates four curves: the gas–water ratio, the contribution of CO₂ to total heat, the contribution of water to total heat, and the cumulative heat extracted, respectively. The local capillary trapping phenomenon can be directly linked to the gas-to-water ratio curve, as shown in Figure 12. Areas with significant local capillary trapping would manifest as stable or show increasing GWR with respect to time, indicating effective CO₂ trapping. Similarly, these trapping zones would impact heat extraction by individual fluids and cumulative heat production, as these regions would contribute less to heat extraction due to trapped CO₂.

4. Discussion

The results depicted in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 furnish a robust, multi-dimensional analysis elucidating the intricate geophysical and thermodynamic interactions underpinning both CO₂ trapping mechanisms and geothermal energy extraction in sedimentary reservoirs. Notably, the exploration into the reservoir heterogeneity is vividly and analytically portrayed through Figure 4 and Figure 8. The detailed exposition of permeability disparities and correlations in these figures is fundamental in predicting the behavior of injected CO₂ [34,37]. It is conspicuous that the variations in permeability and capillary pressures act not merely as passive attributes of the reservoir but as active determinants influencing the mechanics of CO₂ trapping [37,65]. Local capillary trapping, as accentuated, emerges as a notable outcome of these heterogeneities. The effects of such a mechanism are dual-pronged: on the one hand, it catalyzes efficacious CO₂ sequestration, ensuring stable and secure storage [66]; on the other, it poses a potential detriment to heat extraction efficiency, a concern that would be amplified if these trapping zones are in close proximity to thermal recovery sites.

The sequential Figure 5, Figure 6 and Figure 7 serve to chronicle the dynamic responses of the reservoir under disparate operational regimes. Figure 5, delineating the baseline scenario, is juxtaposed against Figure 6, where the injection rate is amplified. This juxtaposition unveils the accelerated migration of supercritical CO₂ seen in the latter scenario, attributable to the elevated flow rates [67]. This brisk advancement of the CO₂ front necessitates an in-depth examination of the temperature and pressure profiles within the reservoir, as these altered conditions bear direct ramifications for the reservoir’s thermodynamic equilibrium [68]. The narrative pivots when Figure 7 enters the discourse, illustrating the scenario wherein the mean permeability is halved. The resulting images display a discernible constriction in the spread of CO₂, alongside moderated temperatures and pressure alterations within the reservoir. This constrained migration and attenuated thermal response underline the profound influence of permeability on reservoir behavior. In these figures, the multifaceted nature of CO₂ storage and geothermal energy extraction is unveiled, elucidating the delicate balance among operational strategies, geological realities, and thermodynamic principles. As each figure appends a layer to our understanding, it becomes unequivocally clear that the success of subsurface interventions is heavily contingent upon a nuanced appreciation of the subsurface’s heterogeneity, coupled with an astute application of geophysical principles and insights.

Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 provide quantifiable metrics—the gas–water ratio (GWR) and cumulative heat extracted (CHE)—which encapsulate the geophysical and thermodynamic complexities into measurable outcomes. Different heterogeneous realizations (R1, R2, etc.) and correlation lengths (25 m, 35 m, and 45 m) manifest as unique GWR curves, confirming that geological heterogeneity directly impacts operational efficiency. The cumulative heat extracted is another metric that provides an integrated view of how much thermal energy is harvested over time, accounting for contributions from both water and scCO₂. The findings align well with the study’s objective, which emphasizes the need for integrated approaches. The quantitative understanding fostered by this work not only strengthens the foundational knowledge for optimizing geothermal and CO₂ storage systems but also provides empirical data that can be leveraged for cross-disciplinary solutions. For instance, these insights into local capillary trapping could be valuable for enhanced oil recovery (EOR) methods, where similar trapping mechanisms are at play.

5. Conclusions

This research was conceived to fulfill the urgent demand for a refined, quantitative comprehension of the convoluted phenomena implicated in CO₂ sequestration and geothermal energy procurement in sedimentary basins. The reservoir parameters are instrumental in dictating the nuances of local capillary trapping, which, in turn, critically affects the efficacy of both CO₂ containment and thermal energy retrieval. Our findings delineate the explicit and measurable consequences of modifying operational parameters such as injection rates and mean permeability for the migration patterns of supercritical CO₂, as well as for the thermal and pressure gradients within the reservoir matrix. Through a longitudinal lens, this inquiry has charted the evolution of state variables over an extended temporal horizon, shedding light on the enduring operational feasibility and the efficiency of energy harnessing and carbon dioxide immobilization strategies. Introducing the gas–water ratio (GWR) and cumulative heat extracted (CHE) as tangible indicators, this study presents a synthesized perspective of the intricate synergies that intertwine geological features, thermodynamic properties, and procedural variables. The ramifications of this research are far-reaching, extending well beyond the immediate sphere of CO₂ storage and geothermics. The deepened understanding of local capillary trapping dynamics and the quantitative indices developed herein could be leveraged to enhance methodologies in cognate disciplines, such as enhanced oil retrieval and the science of subsurface fluid dynamics. While this endeavor has significantly augmented our comprehension, it simultaneously paves the way for subsequent inquiry. Future investigations may profitably delve into additional modalities of geological complexity, examine a spectrum of operational scenarios, and assess the repercussions of geochemical interactions within subsurface conditions.