Impact of Well Placement in the Fractured Geothermal Reservoirs Based on Available Discrete Fractured System

Abstract

Well placement in a given geological setting for a fractured geothermal reservoir is necessary for enhanced geothermal operations. High computational cost associated with the framework of fully coupled thermo-hydraulic-mechanical (THM) processes in a fractured reservoir simulation makes the well positioning a missing point in developing a field-scale investigation. To enhance the knowledge of well placement for different working fluids, we present the importance of this topic by examining different injection-production well (doublet) positions in a given fracture network using coupled THM numerical simulations. Results of this study are examined through the thermal breakthrough time, mass flux, and the energy extraction potential to assess the impact of well position in a two-dimensional reservoir framework. Almost ten times the difference between the final amount of heat extraction is observed for different well positions but with the same well spacing and geological characteristics. Furthermore, the stress field is a strong function of well position that is important concerning the possibility of high-stress development. The objective of this work is to exemplify the importance of fracture connectivity and density near the wellbores, and from the simulated cases, it is sufficient to understand this for both the working fluids. Based on the result, the production well position search in the future will be reduced to the high-density fracture area, and it will make the optimization process according to the THM mechanism computationally efficient and economical.

1. Introduction

Geothermal field development and management is a complex process. Engineering a geothermal system requires appropriate well placement and fracture connectivity to ensure well connectivity and least fluid loss [1,2]. Placement of injection and production wells or a doublet system in a given geological framework to achieve maximum geothermal energy extraction is one of the most complicated and expensive procedures. The location of the injection well concerning production well decides the production mass flux [3,4]. Practically, there is an infinite number of sites where an injection well can be placed in designing an enhanced geothermal system (EGS). Well placement in association with fracture network requires two critical aspects to ensure high heat extraction potential. First, the fractures must be connected sufficiently, and they must provide a high fluid flow rate at a low-pressure difference, and secondly, fluid residence time in the fractures should be increased to allow sufficient heat exchange. Longer residence time enhances the heat extraction capacity and reduces the chances of short-circuiting [5,6]. Figure 1 shows a subsurface fracture network where red is the high-temperature region. Hot water is produced through the red color well from the reservoir, and after passing it through the heat exchanger, it is reinjected to the reservoir with the blue color well. Fractures are the main paths for fluid flow that allow for heat extraction from the various MEET geothermal sites, including Soultz sous Forêts, United Down, Göttingen, and Havelange. Discrete Fracture Network (DFN) characterization is an essential step toward the simulation of reservoir performance. However, the total number of fractures resulting from DFN characterization is a high number (in the order of millions of fractures) which is not feasible for performing numerical simulation (due to computational costs) while considering all of them discretely through the thermal-hydraulic (TH) or thermal-hydraulic-mechanical (THM) simulator. Recently, Lepillier et al. [7] combined TH behavior with the steady-state solid mechanics process to examine the well positioning impact with four doublet scenarios. However, transient temperature and pressure changes will affect the stress field, and four scenarios are not sufficient to accurately demonstrate the impact of well position.

At the same time, the fracture network alignment also contributes to the thermal drawdown, mass flux, and extracted energy. Therefore, it is essential to estimate the well locations a priori for better connectivity and maximum energy extraction. For example, a second well was designed at Rittershoffen (Upper Rhine Graben, France) in the damage zone of the Rittershoffen fault after the drilling of the first one and an additional geophysical survey [8].

This paper considers a two-dimensional fractured reservoir for a potential enhanced geothermal system. Fully coupled thermo-hydro-mechanical (THM) processes are simulated on the fractured reservoir to estimate maximum geothermal energy extraction potential.

Several optimization techniques are available for determining well placement in a reservoir [9]. Some of these methods are gradient-free methods, including genetic algorithms [10], particle swarm optimization algorithm [11], fast marching method [12], and simultaneous perturbation stochastic approximation [13,14], and gradient-based optimization methods, including adjoint methods [15,16,17]. These models lack geological uncertainty, e.g., fracture network connectivity, while considering well placement optimization [18]. Few thermo-hydraulic compositional reservoir simulation-based models on well spacing optimization [19,20,21,22,23,24,25,26]. Based on a coupled thermo-hydraulic model, Akin et al. [27] developed artificial neural networks (ANN) and a search algorithm to optimize an injection well for a geothermal reservoir whereas, Samin et al. [28] developed a hybrid approach integrating a multi-objective genetic algorithm with finite element modeling of thermo-hydraulic processes. EGS involves complex THM processes. Gudmundsdottir and Horne [29] developed an ANN model to characterize fractured geothermal reservoirs for a coupled TH model. Training data necessary for creating a robust ANN model based on a coupled thermo-hydro-mechanical process requires many numerical reservoir simulations. A recent study using an ANN model for a coupled thermo-hydraulic approach supports this idea for a fracture in a hot geothermal reservoir [30] and supercritical geothermal reservoirs [31].

While performing a parametric investigation for a multi-well reservoir, Chen and Jiang [19] found that production well configuration concerning injection well affects the heat mining potential. Chen et al. [32] used a multivariate adaptive regression spline technique coupled with hydrothermal numerical simulation to optimize the well placement under the given fault size and permeability for a prospective geothermal site near Superstition Mountain in Southern California USA. They found that for the maximum net profit over fifty years, the optimal well spacing is 473 m at 30.7 kg/s. For 45° angle between fracture orientation and inlet-outlet connection in a given fracture network, Zhang et al. [33] observed optimized geothermal energy extraction performance. They obtained a stable heat mining rate at reduced efficiency for higher orientation angles. Zhang et al. [2] found that the presence of many fractures in the vicinity of the production well increases the working fluid residence time, and heat recovery efficiency significantly improves. They suggested that thermally-induced fractures near the production well assist in greater power generation than when the fracture density is high in the vicinity of the injection well. They also observed that placing the production well in the high permeability region increases heat production. Gao et al. [34] used a coupled thermo-hydraulic model for a discrete fracture network in a fractured geothermal reservoir to investigate heat extraction performance. They used multilateral well orientations with a varying number of branch wells and well orientation. They found that production temperature decreases with an increase in the well and fracture intersections, whereas injection pressure increases. Aliyu et al. [35] and Aliyu and Chen [36] used COMSOL Multiphysics to develop a model depicting THM and TH processes in a geothermal reservoir for two fractures and single fractures, respectively. They estimated the impact of well spacing on thermal energy extraction performance.

The MEET project framework considers water as the working fluid or heat-carrying fluid from a geothermal system. However, this study finds CO₂ as an alternative to water because the loss of CO₂ as the heat-carrying fluid is environment friendly [37,38]. Furthermore, the use of supercritical CO₂ may assist in the formation of an interconnected fracture network of multiple channels at a lower pressure than water due to smaller fluid density and viscosity [38]. Due to the lower reactivity of CO₂ in comparison to water, the possible silica dissolution and precipitation at high temperatures and pressure decreases [37,38,39]. Additionally, Bongole et al. [40] observed that the reservoir deformation is more minor when CO₂ is the working fluid compared to water due to the lower heat capacity of CO₂. The lower freezing point of CO₂ than water helps heat rejection at a much lower temperature. Therefore, its geothermal systems may work even for cold climatic conditions where water is unusable [41]. The above literature shows no available THM model for determining the well placement in a geothermal reservoir concerning a given fracture map to maximize mass flux, energy extraction, and the thermal drawdown duration. In this work, a fully coupled THM model is developed and used to demonstrate the importance of well position by characterizing the fracture network connectivity and density. This study will build a basis for future well placement optimization considering THM processes for a given fracture network. The present study is organized in the following manner. First, a mathematical and numerical model is presented for a coupled THM process followed by results and discussion on optimizing well positions in a two-dimensional fracture network based on thermal drawdown, mass flux, and energy extraction potential followed by conclusions.

2. Methodology

This study uses a fracture network based on outcrop fractures mapped from Otsego County in New York state [42] for THM modeling. The total number of fractures in this outcrop map is 440. As depicted in Figure 1, this study considers a two-dimensional geometry at subsurface conditions. The reservoir geometry is a two-dimensional planar model (1000 m × 600 m), and the injection-production wells are placed 500 m apart. An initial case is considered where the injection well is present at point 0 and the production well at point 180, as shown in Figure 2. Considering this axis as diameter, a circular zone is assumed, and the perimeter is divided into 36 equal intervals. These 36 intervals are considered for placing the injection and production wells, and they are arranged at α angle from the base case. All fractures are assumed as interior boundaries, and the displacement is constrained in all normal directions. The side boundaries are assumed as no flow for both heat and mass exchange.

Conservation equation of mass when coupled with pore volume and fluid temperature alteration for a porous medium is [43]:

$${\rho }_1\left({\varphi }_m{S}_1+\left(1-{\varphi }_m\right)S_m\right)\frac{\partial p}{\partial t}-{\rho }_1\left({\alpha }_m\left({\varphi }_m{\beta }_1+\left(1-{\varphi }_m\right){\beta }_m\right)\right)\frac{\partial T}{\partial t}+{\rho }_1{\alpha }_m\frac{\partial {\epsilon }_V}{\partial t}=\nabla ·\left(\frac{{\rho }_1{k}_m}{\mu }\nabla p\right)$$

(1)

All the parameters are listed in Appendix A. Water and supercritical CO₂ are considered heat-transmitting fluids in this study. The equation that governs fluid flow along the internal fractures is:

$${\rho }_1\left({\varphi }_f{S}_1+\left(1-{\varphi }_f\right)S_{mf}\right)e_h\frac{\partial p}{\partial t}-{\rho }_1\left({\alpha }_f\left({\varphi }_f{\beta }_1+\left(1-{\varphi }_f\right){\beta }_f\right)\right)e_h\frac{\partial T}{\partial t}+{\rho }_1{\alpha }_f{e}_h\frac{\partial {\epsilon }_V}{\partial t}={\nabla }_T·\left(\frac{e_h{\rho }_1{k}_f}{\mu }{\nabla }_{T}p\right)+n·Q_m$$

(2)

In Equation (2), fluid flow along the fracture width is ignored because fracture aperture is much smaller than the fracture length. Fractures and the rock matrix are assumed at thermodynamic inequilibrium. In other words, the local thermal non-equilibrium model is implemented in this investigation.

$$\left(1-{\varphi }_m\right){\rho }_m{C}_{p,m}\frac{\partial T_m}{\partial t}=\nabla ·\left(\left(1-{\varphi }_m\right){\lambda }_m\nabla T_m\right)+q_{ml}\left(T_l-T_m\right)$$

(3)

$$\left(1-{\varphi }_f\right)e_h{\rho }_f{C}_{p,f}\frac{\partial T_m}{\partial t}={\nabla }_T·\left(\left(1-{\varphi }_f\right)e_h{\lambda }_f{\nabla }_T{T}_m\right)+e_h{q}_{fl}\left(T_l-T_m\right)+n·\left(-\left(1-{\varphi }_m\right){\lambda }_m\nabla T_m\right)$$

(4)

The energy balance equation for the rock matrix and fractures are shown by Equations (3) and (4), respectively. The energy balance equation for either water or CO₂ is:

$${\varphi }_m{\rho }_l{C}_{p,l}\frac{\partial T_l}{\partial t}+{\varphi }_m{\rho }_l{C}_{p,l}\left(-\frac{k_m\nabla p}{\mu }\right)·\nabla T_l=\nabla .\left({\varphi }_m{\lambda }_l\nabla T_l\right)+q_{ml}\left(T_m-T_l\right)$$

(5)

The following equation can write heat exchange between a rock and the fracture matrix:

$${\varphi }_f{e}_h{\rho }_l{C}_{p,l}\frac{\partial T_l}{\partial t}+{\varphi }_f{e}_h{\rho }_l{C}_{p,l}\left(-\frac{k_f{\nabla }_{T}p}{\mu }\right)·{\nabla }_T{T}_l={\nabla }_T·\left({\varphi }_f{e}_h{\lambda }_l{\nabla }_T{T}_l\right)+e_h{q}_{fl}\left(T_m-T_l\right)+n·\left(-{\varphi }_l{\lambda }_l\nabla T_l\right)$$

(6)

In Equation (6), the Darcy flux in the fractures is $u_f=-\frac{k_f{\nabla }_{T}p}{\mu }$ and heat flux is $n.q_l=n·\left(-{\varphi }_l{\lambda }_l\nabla T_l\right)$.

A fully coupled thermo-hydro-mechanical model is developed in this study. If effective stress is ${\sigma }_{eff}^{ij}={\sigma }_{ij}+{\alpha }_{p}p{\delta }_{ij}$ and the volumetric expansion coefficient of porous media is ${\beta }_T={\varphi }_l{\beta }_l+\left(1-{\varphi }_m\right){\beta }_m$, then the stress-strain relationship considering fully coupled thermoelastic and poroelastic stress can be written as:

$${\sigma }_{ij}=2G{\epsilon }_{ij}+\lambda tr\epsilon {\delta }_{ij}-{\alpha }_{p}p{\delta }_{ij}-K^{\prime }{\beta }_{T}T{\delta }_{ij}$$

(7)

The reservoir deformation equation can be written as:

$$G{u}_{i,jj}+\left(G+\lambda \right)u_{j,ji}-{\alpha }_p{p}_{,i}-K^{\prime }{\beta }_T{T}_{,i}+f_i=0$$

(8)

The opening and closure of the thermo-poroelastic stress-dependent fracture aperture are modeled using the Barton and Bandis model [44,45] as follow:

$$\Delta e_n=\frac{e_0}{1+9\frac{{\sigma }_{eff}^n}{{\sigma }_{nref}}}$$

(9)

In Equation (9), $\Delta e_n$ is the fracture aperture change under in-situ stress conditions.

Thermodynamic properties of water and CO₂ are represented by dynamic viscosity (Equations (10) and (11)), specific heat capacity (Equations (12) and (13)), density (Equations (14) and (15)), and thermal conductivity (Equations (16) and (17)) [42], and they are implemented in Equations (1)–(6).

$${\mu }_w=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 K\right)$$

(10)

$${\mu }_{C{O}_2}=-1.49\times {10}^{-6}-6.47\times {10}^{-8}\times T^1-3.66\times {10}^{-11}\times T^2+1.25\times {10}^{-14}\times T^3 \left(220-600 K\right)$$

(11)

$$C_{p,w}=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$$

(12)

$$C_{p,C{O}_2}=459.91+1.86\times T^1-2.13\times {10}^{-3}\times T^2+1.22\times {10}^{-6}\times T^3 \left(220-600 K\right)$$

(13)

$${\rho }_w=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)

$${\rho }_{C{O}_2}=pA\times 0.04401/RT$$

(15)

$${\kappa }_w=-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$$

(16)

$${\kappa }_{C{O}_2}=-1.32\times {10}^{-3}+4.14\times {10}^{-5}\times T^1+6.71\times {10}^{-8}\times T^2-2.11\times {10}^{-11}\times T^3$$

(17)

In Equation (15), $pA$ is the absolute pressure, and $R$ is the molar gas constant. Coefficients in Equations (10)–(17) are constants and obtained from various correlations [42].

COMSOL Multiphysics version 5.5 [43] is used to perform numerical modeling of THM processes. It uses a finite element method to solve general purpose partial differential equations. The full mesh contains 112,818 domain elements and 13,071 boundary elements. This free triangular mesh is generated by using the maximum element size of 37 m, and minimum element size of 0.125 m, maximum element growth of 1.25 with the curvature factor 0.25, and the resolution of the narrow regions is 1. For the numerical modeling purpose, we have used a scaled absolute tolerance of magnitude 10⁻⁸ and automatic time step constraint. We assumed Backward Differentiation Formula (BDF) for time-stepping with maximum BDF order as two and minimum BDF order as 1. Further, we have validated our model with a soil thermal consolidation model as Bai [46] demonstrated in Mahmoodpour et al. [47].

3. Results and Discussions

Numerical simulation results from coupled THM mechanisms associated with a geothermal energy extraction process from a fractured reservoir are presented in this section. First, we performed a sensitivity analysis for three different mesh elements. For the following stage, the adopted sequence of presentation is:

(a)

coupled THM mechanisms for heat mining using water as heat-carrying fluid,

(b)

coupled THM processes when CO₂ is the heat-carrying fluid, and

(c)

predicting a suitable doublet well position for a given fracture network to obtain highest mass flux from the production well and maximize the heat production.

The results are presented for two working fluids: water and CO₂. Reservoir permeability of 2 mD and 5 mD are considered. These values are chosen in a way that sweep efficiency with the different working fluids to be similar at the same time.

Furthermore, permeability values are kept constant to understand the working fluid effect by running another case 180. Other parameters are listed in

Table 1. Numerical simulation parameters (see the range of database in [48]).

Parameter	Magnitude for Water-Based Simulations	Magnitude for CO2 Based Simulations
Young’s modulus	$40 \mathrm{GPa}$	$40 \mathrm{GPa}$
Poisson’s ratio	$0.25$	$0.25$
Rock density	$2500 \frac{\mathrm{kg}}{{\mathrm{m}}^3}$	$2500 \frac{\mathrm{kg}}{{\mathrm{m}}^3}$
Horizontal stress	$50 \mathrm{MPa}$	$50 \mathrm{MPa}$
Vertical stress	$50 \mathrm{MPa}$	$50 \mathrm{MPa}$
Initial pressure	$30 \mathrm{MPa}$	$30 \mathrm{MPa}$
Injection pressure	$50 \mathrm{MPa}$	$50 \mathrm{MPa}$
Rock porosity	$0.2$	$0.2$
Rock permeability	$2 \mathrm{mD}$	$5 \mathrm{mD}$
Fracture zone porosity	$0.5$	$0.5$
Fracture roughness	$1$	$1$
Fracture aperture	$0.2 \mathrm{mm}$	$0.2 \mathrm{mm}$
Closure stress	$150 \mathrm{MPa}$	$150 \mathrm{MPa}$
Wellbore radius	$0.2 \mathrm{m}$	$0.2 \mathrm{m}$
Rock thermal conductivity	$3 \frac{\mathrm{W}}{\mathrm{m}\times \mathrm{K}}$	$3 \frac{\mathrm{W}}{\mathrm{m}\times \mathrm{K}}$
Fracture zone thermal conductivity	$2.5\frac{\mathrm{W}}{\mathrm{m}\times \mathrm{K}}$	$2.5\frac{\mathrm{W}}{\mathrm{m}\times \mathrm{K}}$
Rock specific heat capacity	$800\frac{\mathrm{J}}{\mathrm{kg}\times \mathrm{K}}$	$800\frac{\mathrm{J}}{\mathrm{kg}\times \mathrm{K}}$
Fracture zone specific heat capacity	$800\frac{\mathrm{J}}{\mathrm{kg}\times \mathrm{K}}$	$800\frac{\mathrm{J}}{\mathrm{kg}\times \mathrm{K}}$
Initial temperature	$200 °\mathrm{C}$	$200 °\mathrm{C}$
Biot coefficient	$0.7$	$0.7$
Thermal expansion coefficient	${10}^{-5}\frac{1}{\mathrm{K}}$	${10}^{-5}\frac{1}{\mathrm{K}}$
Injection temperature	$70 °\mathrm{C}$	$70 °\mathrm{C}$

that are not site specific and selected to represent a generic geothermal system, and the fracture map is the same for all scenarios. Constant injection and production pressures are considered at both the wellbores. A two-dimensional horizontal cross-sectional reservoir is considered for all the simulations.

We performed a mesh sensitivity analysis with case 180 for CO₂ (it has the highest velocity variation, and the simulations convergence is the most complex among all the cases). Here, the mesh sensitivity analysis is attached for the simulations with (a) 92,655 domain elements and 12,103 boundary elements, (b) 112,818 domain elements and 13,071 boundary elements and (c) 181,410 domain elements and 15,687 boundary elements. Convergence was not achieved with the mesh size of the 71,089 domain elements and 10,357 boundary elements. The maximum element size for the standard case is (a) 67 m, (b) 37 m, and (c) 18 m, whereas the minimum element size is (a) 0.3 m, (b) 0.125 m, and (c) 0.075 m. The maximum element growth rate is (a) 1.3, (b) 1.25, and (c) 1.2, the curvature factor is (a) 0.3, (b) 0.25, and (c) 0.25, and the resolution of the narrow regions is 1 for all three cases. Free triangular meshes are used for discretizing this domain.

Based on this description, our results are insensitive to the mesh refinements (see Figure 3).

Figure 4 shows the time evolution of reservoir temperature distribution during the heat extraction operation using water and CO₂ for case 180. Here, the injection well is present at 180, and the production well is present at 0 as shown in Figure 2. The reservoir permeability for the left (Figure 4(a1–d1)) and right (Figure 4(a3–d3)) columns is 5 mD whereas the middle column has 2 mD permeability (Figure 4(a2–d2)). Higher reservoir permeability in the case of the left column causes faster cold fluid propagation through the fracture network. Additionally, water propagation through the fractures becomes less dominant, and it starts flowing through the rock matrix at higher permeability as shown in Figure 4(c1,d1). We adopted smaller reservoir permeability for well placement when water is the working fluid to account for both these factors. The cold-water propagation is aligned along the dominant fracture rather than the horizontal axis between the doublet. The reason behind selecting different permeability values for water and CO₂ is to reach a similar sweep efficiency with different fluids; however, we provided the quantitative comparison between simulation of water at 5 mD, 2 mD and CO₂ at 5 mD for case 180 as shown in Figure 5. Figure 5a–c shows that CO₂ is a better-working fluid concerning the breakthrough time, mass flux, and cumulatively extracted energy, respectively over 30 years.

The viscosity of supercritical CO₂ at injection conditions is approximately half compared to water. Higher reservoir permeability and lower viscosity indicate CO₂ propagation through the fractures as well as through the matrix rather than flowing through fractures only as seen in Figure 4(a2,b2). Figure 4(a3–d3) shows that the cold fluid plume spread is much diffusive compared to water, and flow is primarily occurring through the matrix. However, Figure 4(d3) shows that CO₂ flows through the dominant fractures near the production well. Therefore, fluid propagation through the fractures is the principal mechanism between the doublet, which is assisted by flow through permeable rock matrix.

Furthermore, convective heat transfer inside the rock matrix and the fracture is the primary heat transfer mechanism. Therefore, the fractures control the heat transfer in a fractured reservoir. To show the relative importance of convective to conductive heat transfer, we calculated Peclet number ($Pe$), a nondimensional number which indicates the ratio of convective to conductive heat transfer [49] and it can be written as $Pe=\frac{uL{\rho }_i{C}_{p,i}}{k_i}$, where $u$ is the fluid velocity, $L$ is the characteristic length (here it is 500 m, the distance between the two wells), the subscript $i$ indicates the fluid either water or CO₂, ${\rho }_i$ is the fluid density, $C_{p,i}$ is the specific heat capacity of fluid and $k_i$ is the fluid thermal conductivity. Figure 6a shows the $Pe$ value after one year for the entire reservoir for case 180, where CO₂ is the working fluid. To elaborate on the relative impact of convective heat transfer inside the fracture, one fracture is selected as shown in Figure 6b, and the corresponding $Pe$ number is shown in Figure 6c. $Pe$ number is estimated for five different times, and for all the cases, $Pe$ number is significantly larger than 1, indicating more vigorous convection than conduction.

Results obtained from the 36 reservoir simulations on the well positioning are shown in Figure 7. For water simulations, reservoir permeability is 2 mD, whereas, for cases with CO₂ as working fluid, permeability is 5 mD. Water-based models demonstrate faster thermal breakthrough (Figure 7a) due to higher specific heat capacity than CO₂. Therefore, water simulations are presented for 30 years whereas CO₂ results are plotted for 300 years. We magnified the results for the CO₂ over 30 years in Figure 7b,d,f to compare it with water (see Figure 7a,c,e). Figure 7a shows thermal drawdown at the production well when water is the working fluid. The fastest thermal drawdown was observed for case 220 (see Figure 8(c1)) whereas the slowest thermal drawdown occurred for case 130 (see Figure 8(a1)). From Figure 2, it is clear that case 220 has a well position along a dominant fracture supported by minor intersecting fractures, whereas case 130 wells are aligned approximately orthogonal to this prevalent fracture. Figure 8(a2,a3) show thermoelastic stress along the horizontal and vertical directions, respectively, indicating stress localization spans across the cold fluid plume region. Greater concentration of connected fractures in the area away from the doublet axis causes prolonged thermal breakthrough time. Therefore, in 30 years, the temperature drop is approximately 40 °C for case 130, whereas case 220 shows a 75 °C temperature drop at the production well.

In comparison to water, CO₂ has approximately seven times smaller thermal conductivity at the injection conditions. Due to this, thermal depletion time is prolonged when CO₂ is the working fluid compared to water as the heat-carrying fluid. Figure 7b shows the thermal drawdown at the production well when CO₂ is the operating fluid. It shows that the thermal drawdown curves depend significantly on the fracture network connectivity than water. The slowest thermal drawdown is demonstrated by case 40, where production well temperature drops by approximately 20 °C in 300 years. In contrast, the fastest thermal drawdown is displayed by case 180, where around 90 °C temperature drop is estimated. Figure 9(a1,a2) show reservoir temperature distribution for case 40 and case 180, respectively. In Figure 9(a1), the cold fluid spread is extremely slow since a high fracture density is present near point 40, as shown in Figure 2 that is present away from the doublet axis. This leads to a reduced amount of cold fluid injection and restricted heat exchange between the fluid–fractures and fluid–matrix in the reservoir, decreasing the horizontal and vertical thermoelastic stress as shown in Figure 9(a2,a3), respectively. A detailed sensitivity analysis of dependent parameters is performed for water [47] and CO₂ [50] based geothermal systems for the same fractured reservoir as mentioned in this paper. Our numerical simulations consider poroelastic stress, but we have not shown here that they contribute little due to the fluid injection and production, as shown in our previous findings [47,50]. For case 180, Figure 9(b1) shows reservoir temperature distribution after 300 years. It indicates that the hot fluid has been completely extracted between the doublet, and the heat replenishment is too slow to recharge this depleting heat content. Figure 9(b2,b3) approves this reasoning that due to favorable fracture density along the doublet axis, higher fluid flux reinjection results in higher thermoelastic stress evolution.

Figure 7c,d show mass flux at the production well for 36 cases when the working fluid in the reservoir is water and CO₂, respectively. The mass flux for CO₂ is approximately five times higher than water to compensate for smaller viscosity and higher permeability by maintaining the reservoir injection pressure. For both the fluids highest mass flux is observed for case 180, and the smallest mass flux is marked for case 350. Even though these two doublet arrangements have approximately the same axis (endpoints of a single line connecting injection and production wells), the fracture density near the production well plays a crucial role in mass flux; a similar observation was made by Zhang et al. [2]. For case 180, fractures are well connected near the production well, which assists in higher fluid production, whereas in the case of 350, fractures are not connected in a wide area leading to smaller fluid production. The temperature front in Figure 8(e1) shows the weak convective flow for case 350. This can be easily seen from the stress distribution plots in Figure 8(e2,e3) for case 350 when water is the working fluid and in Figure 9(c2,c3) for case 350 when CO₂ is the working fluid. The decrease in mass flux for all the cases with time is due to a in water viscosity with increased fluid temperature. However, we observe that the mass flux increases with time if CO₂ is the working fluid. This increase is approximately <30% between the period when CO₂ production starts till 300 years of numerical simulation. This increase is pronounced for the case 180 where we observe that the mass flux increases from 1.15 to 1.5 kg/s and the increase starts after approximately 50 years from the beginning of the operation. This discrepancy is observed due to limitations in the equation of state used in modeling using COMSOL Multiphysics. Since viscosity is a function of temperature only, the mass flux increase is observed after the breakthrough time (see Equation (11)).

The energy extraction potential from the reservoir for both the fluids are approximately the same (see Figure 7e,f) since in the case of water simulations, reservoir permeability is 2.5 times more negligible compared to CO_2, and the mass flux of CO₂ is approximately 3–5 times greater than water. In contrast, the specific heat capacity of water is around three times higher compared to CO₂ with three times higher viscosity. Due to higher mass flux and delayed thermal drawdown operation, total energy extraction potential is significantly higher when CO₂ is the working fluid. Case 180 for both the fluids shows the highest energy extraction potential due to the maximum mass flow rate in Figure 8(b1–b3). However, due to higher mass flux, thermal depletion is also fastest, and therefore, doublet placed for case 180 may not show a longer operation when water is the working fluid. Case 250 for water and case 350 for CO₂ show the least energy extraction potential over 30 and 300 years, respectively. Figure 8(d1) shows the reservoir temperature distribution for case 250 when water is used for heat transmission. It is visible from Figure 8d1 that a more incredible amount of cold fluid is present near the injection well, and there is only one large fracture along the doublet axis. This limits the fluid transmission at a higher rate only through a narrow region causing limited energy extraction.

Furthermore, Figure 8(d2,d3) shows the corresponding stress distribution for the horizontal and vertical directions. On the other hand, Figure 9(c1) shows reservoir temperature distribution for case 350 and CO₂ is used as working fluid and localization of cold fluid near the injection well in the absence of any dominating fracture system. The passage of fluid is limited through the fractures toward the production well. The stress field in Figure 8(c2,c3) shows the horizontal and vertical stress are well aligned with the temperature propagation. However, since thermal breakthrough is slower when CO₂ is used for heat transmission, energy extraction potential may enhance if EGS operation is performed beyond 300 years.

4. Conclusions

Geothermal energy extraction from deep fractured reservoirs can support high energy demand for a long duration. Water and CO₂ are two fluids that can extract energy from the subsurface. A fractured reservoir shows a complex network of fractures, and fracture conductivity controls the primary fluid passage for heat extraction longevity of the operation. Well placement for a given fracture network should consider the fracture density and orientation. Keeping all the parameters constant except the injection–production doublet axis orientation, we observe a difference of approximately ten times of energy extracted among the studied cases. High fracture density in the vicinity of the production well is the reason behind this increased energy extraction. The doublet axis orientation affects the injectivity (poroelastic stress) and temperature propagation (thermoelastic stress). It has a great impact on the stress field development during heat extraction.

Fluid type plays a significant role in determining the THM behavior of the EGS. The viscosity of fluid determines the temperature propagation through the fractures, as well as through the rock matrix. CO₂ with lower viscosity can penetrate easily inside the matrix zone. This effect, combined with the lower specific heat capacity of CO₂, eventuates the cold front of fluid propagation through the matrix and fracture. While water with high viscosity and specific heat capacity mainly transmits heat alongside the fracture and results in early breakthrough time. Different cases with water have a small range of breakthrough time compared to CO₂. While CO₂ shows a higher flow rate, resulting from the lower viscosity, this behavior is compensated by the higher heat capacity of water. Therefore, the overall heat extraction is comparable for both fluids.