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

: 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 ﬁeld-scale investigation. To enhance the knowledge of well placement for different working ﬂuids, 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 ﬂux, 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 ﬁnal amount of heat extraction is observed for different well positions but with the same well spacing and geological characteristics. Furthermore, the stress ﬁeld 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 sufﬁcient to understand this for both the working ﬂuids. 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 efﬁcient and economical.


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 lowpressure 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.
Geosciences 2021, 11, x FOR PEER REVIEW 2 of 19 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 steadystate 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.
Figure 1.An enhanced geothermal system.The subsurface plan shows an intricate network of fractures.For optimized power generation, appropriate placement of the injection and production wells is necessary.
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 Production well Injection well 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 2 as an alternative to water because the loss of CO 2 as the heat-carrying fluid is environment friendly [37,38].Furthermore, the use of supercritical CO 2 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 2 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 2 is the working fluid compared to water due to the lower heat capacity of CO 2 .The lower freezing point of CO 2 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.

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.

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]: All the parameters are listed in Appendix A. Water and supercritical CO2 are considered heat-transmitting fluids in this study.The equation that governs fluid flow along the internal fractures is: 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.Conservation equation of mass when coupled with pore volume and fluid temperature alteration for a porous medium is [43]: All the parameters are listed in Appendix A. Water and supercritical CO 2 are considered heat-transmitting fluids in this study.The equation that governs fluid flow along the internal fractures is: 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. (1 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 2 is: The following equation can write heat exchange between a rock and the fracture matrix: In Equation ( 6), the Darcy flux in the fractures is A fully coupled thermo-hydro-mechanical model is developed in this study.If effective stress is σ ij e f f = σ ij + α p pδ ij and the volumetric expansion coefficient of porous media is then the stress-strain relationship considering fully coupled thermoelastic and poroelastic stress can be written as: The reservoir deformation equation can be written as: The opening and closure of the thermo-poroelastic stress-dependent fracture aperture are modeled using the Barton and Bandis model [44,45] as follow: In Equation ( 9), ∆e n is the fracture aperture change under in-situ stress conditions.Thermodynamic properties of water and CO 2 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).
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 −8 and automatic time step constraint.We assumed Backward Differentiation Formula (BDF) for timestepping 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].

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 2 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 2 .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 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 2 (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 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 2 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 2 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 2 at 5 mD for case 180 as shown in Figure 5. Figure 5a-c shows that CO 2 is a better-working fluid concerning the breakthrough time, mass flux, and cumulatively extracted energy, respectively over 30 years.

Injection temperature 70 °C 70
We performed a mesh sensitivity analysis with case 180 for CO2 (it velocity variation, and the simulations convergence is the most complex cases).Here, the mesh sensitivity analysis is attached for the simulations domain Based on this description, our results are insensitive to the mesh r Figure 3).Figure 4 shows the time evolution of reservoir temperature distribu heat extraction operation using water and CO2 for case 180.Here, the i The reason behind selecting different permeability values for water and CO2 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 CO2 at 5 mD for case 180 as shown in Figure 5. Figure 5a-c shows that CO2 is a better-working fluid concerning the breakthrough time, mass flux, and cumulatively extracted energy, respectively over 30 years.The viscosity of supercritical CO2 at injection conditions is approximately half compared to water.Higher reservoir permeability and lower viscosity indicate CO2 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 The viscosity of supercritical CO 2 at injection conditions is approximately half compared to water.Higher reservoir permeability and lower viscosity indicate CO 2 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 2 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 = , 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 2 , ρ i is the fluid density, C p,i is the 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 2 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 2 as working fluid, permeability is 5 mD.Water-based models demonstrate faster thermal breakthrough (Figure 7a) due to higher specific heat capacity than CO 2 .Therefore, water simulations are presented for 30 years whereas CO 2 results are plotted for 300 years.We magnified the results for the CO 2 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 2 has approximately seven times smaller thermal conductivity at the injection conditions.Due to this, thermal depletion time is prolonged when CO 2 is the working fluid compared to water as the heat-carrying fluid.Figure 7b shows the thermal drawdown at the production well when CO 2 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 2 [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.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 CO2 as working fluid, permeability is 5 mD.Water-based models demonstrate faster thermal breakthrough (Figure 7a) due to higher specific heat capacity than CO2.Therefore, water simulations are presented for 30 years whereas CO2 results are plotted for 300 years.Figure 7c,d show mass flux at the production well for 36 cases when the working fluid in the reservoir is water and CO2, respectively.The mass flux for CO2 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 CO2 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 CO2 is the working fluid.This increase is approximately <30% between the period when CO2 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 Figure 7c,d show mass flux at the production well for 36 cases when the working fluid in the reservoir is water and CO 2 , respectively.The mass flux for CO 2 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 2 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 2 is the working fluid.This increase is approximately <30% between the period when CO 2 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 2 is approximately 3-5 times greater than water.In contrast, the specific heat capacity of water is around three times higher compared to CO 2 with three times higher viscosity.Due to higher mass flux and delayed thermal drawdown operation, total energy extraction potential is significantly higher when CO 2 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 2 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 2 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 2 is used for heat transmission, energy extraction potential may enhance if EGS operation is performed beyond 300 years.

Conclusions
Geothermal energy extraction from deep fractured reservoirs can support high energy demand for a long duration.Water and CO 2 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 2 with lower viscosity can penetrate easily inside the matrix zone.This effect, combined with the lower specific heat capacity of CO 2 , 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 2 .While CO 2 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.Technische Universität Darmstadt has provided institutional support to authors.Authors would like to acknowledge anonymous reviewers for their valuable comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.An enhanced geothermal system.The subsurface plan shows an intricate network of fractures.For optimized power generation, appropriate placement of the injection and production wells is necessary.

Figure 2 .
Figure 2. The geometry of the reservoir.Injection and production wells are placed 500 m apart.Total 36 cases or 36 values of α are considered for simulations.Here, when the injection well is present at 0 and production well is present at 180 and they are 500 m placed apart, the value of α is 0. When the injection well is present at 340 and the production well is present at 160, the value of α is 340.

Figure 2 .
Figure 2. The geometry of the reservoir.Injection and production wells are placed 500 m apart.Total 36 cases or 36 values of α are considered for simulations.Here, when the injection well is present at 0 and production well is present at 180 and they are 500 m placed apart, the value of α is 0. When the injection well is present at 340 and the production well is present at 160, the value of α is 340.
(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.
elements and 12,103 boundary elements, (b) 112,818 domain elem boundary elements and (c) 181,410 domain elements and 15,687 boundary vergence was not achieved with the mesh size of the 71,089 domain elem boundary elements.The maximum element size for the standard case is m, and (c) 18 m, whereas the minimum element size is (a) 0.3 m, (b) 0.125 m.The maximum element growth rate is (a) 1.3, (b) 1.25, and (c) 1.2, the c is (a) 0.3, (b) 0.25, and (c) 0.25, and the resolution of the narrow regions cases.Free triangular meshes are used for discretizing this domain.

Figure 3 .
Figure 3. Mesh sensitivity results for three different number of mesh elements.

Figure 3 .
Figure 3. Mesh sensitivity results for three different number of mesh elements.

Figure 4 .
Figure 4. Reservoir temperature distribution at time 1 year (a1-a3), 5 years (b1-b3), 10 years (c1-c3) and 30 years (d1-d3) when the injection well is present at 180 and the production well is placed at 0. Results from water for reservoir permeability 5 mD and 2 mD are shown by (a1-d1) and (a2-d2) respectively.CO2 as working fluid results are displayed in (a3-d3).The reservoir permeability for CO2 simulations is 5 mD.The injection wellbore position is shown by cross symbol whereas a circle indicates production wellbore position.

Figure 4 .Figure 5 .
Figure 4. Reservoir temperature distribution at time 1 year (a1-a3), 5 years (b1-b3), 10 years (c1-c3) and 30 years (d1-d3) when the injection well is present at 180 and the production well is placed at 0. Results from water for reservoir permeability 5 mD and 2 mD are shown by (a1-d1) and (a2-d2) respectively.CO 2 as working fluid results are displayed in (a3-d3).The reservoir permeability for CO 2 simulations is 5 mD.The injection wellbore position is shown by cross symbol whereas a circle indicates production wellbore position.Geosciences 2021, 11, x FOR PEER REVIEW 9 of 19

Figure 5 .
Figure 5.Comparison of working fluid effects, (a) temperature at the production well, (b) mass flux of different fluid, and (c) extracted energy.

Geosciences 2021 , 19 Figure 6 .
Figure 6.For case 180 and CO2 as the working fluid after one year, (a) Peclet number distribution across the reservoir, (b) fracture considered for estimating fracture length, and (c) Peclet number along with the fracture length.

Figure 6 .
Figure 6.For case 180 and CO 2 as the working fluid after one year, (a) Peclet number distribution across the reservoir, (b) fracture considered for estimating fracture length, and (c) Peclet number along with the fracture length.

Figure 7 .
Figure 7.The temperature at the production well for (a) water and (b) CO2, mass flux at the production well for (c) water and (d) CO2, and cumulative energy extraction using (e) water and (f) CO2 as heat-carrying fluid.Results from 36 simulation cases are plotted as shown by the legend, and case 180 is indicated by bold magenta color.Here Case 180 means injection well is present at 180, and production well is present at 0 (see Figure 2).

Figure 7 .
Figure 7.The temperature at the production well for (a) water and (b) CO 2 , mass flux at the production well for (c) water and (d) CO 2 , and cumulative energy extraction using (e) water and (f) CO 2 as heat-carrying fluid.Results from 36 simulation cases are plotted as shown by the legend, and case 180 is indicated by bold magenta color.Here Case 180 means injection well is present at 180, and production well is present at 0 (see Figure 2).
flux exchange between porous media and the fracture ∇ T Gradient operator restricted to the fracture's tangential plane Specific heat capacity of the rock matrix λ m Heat conductivity of the rock matrixTable A1.Cont.Symbol Parameter q ml Rock matrix-pore fluid interface heat transfer coefficient ρ f density of the fracture zone C p, f Specific heat capacity of the fracture λ f Heat conductivity of the fracture q f l Rock fracture-fluid interface heat transfer coefficient C p & C p,l Heat capacity of the fluid at a constant pressure 2ν) , bulk modulus of the drained porous media β T Volumetric thermal expansion coefficient of porous media Change in the initial aperture of the fracture under in-situ stresses Effective normal stress acting on the fracture surface σ nre f Effective normal stress required to cause 90% reduction in fracture aperture µ CO 2 dynamic viscosity κ CO 2 thermal conductivity