*Article* **Numerical Investigation on the Influence of Surface Flow Direction on the Heat Transfer Characteristics in a Granite Single Fracture**

**Xuefeng Gao 1, Yanjun Zhang 1,2, Zhongjun Hu 1,\* and Yibin Huang <sup>1</sup>**


**Abstract:** As fluid passes through the fracture of an enhanced geothermal system, the flow direction exhibits distinct angular relationships with the geometric profile of the rough fracture. This will inevitably affect the heat transfer characteristics in the fracture. Therefore, we established a hydrothermal coupling model to study the influence of the fluid flow direction on the heat transfer characteristics of granite single fractures and the accuracy of the numerical model was verified by experiments. Results demonstrate a strong correlation between the distribution of the local heat transfer coefficient and the fracture morphology. A change in the flow direction is likely to alter the transfer coefficient value and does not affect the distribution characteristics along the flow path. Increasing injection flow rate has an enhanced effect. Although the heat transfer capacity in the fractured increases with the flow rate, a sharp decline in the heat extraction rate and the total heat transfer coefficient is also observed. Furthermore, the model with the smooth fracture surface in the flow direction exhibits a higher heat transfer capacity compared to that of the fracture model with varying roughness. This is attributed to the presence of fluid deflection and dominant channels.

**Keywords:** heat transfer; granite fracture; flow direction; enhanced geothermal system

#### **1. Introduction**

Geothermal energy is a clean and environmentally friendly renewable energy source with a wide distribution range, large reserves, and long duration [1–3]. The high temperature rock mass buried 3–10 km underground is characterized by low permeability and low porosity [4]. In order to extract thermal energy from such high temperature rock mass, several countries (led by the United States) proposed the use of artificial hydraulic fracturing, which led to the development of reservoir-transforming hydro-shearing techniques to establish fracture network channels and improve the heat transfer capacity [5]. This geothermal project is often referred to as an enhanced geothermal system (EGS). The fracture surface formed by hydraulic fracturing is typically rough, which directly affects the heat transfer performance of the working medium (water or CO2) [6].

The key to extracting geothermal resources by EGS is the flow and heat exchange process of working fluid in fractures of high temperature rock mass [7]. However, when the fluid begins to flow from the injection well into the fracture, the initial direction of flow often has an angle equal to the direction of the rough fracture surface, which inevitably affects the heat exchange effect. The accurate determination of the convection heat transfer coefficient between the fluid and fracture is crucial to the optimization of the reservoir transformation and productivity predictions [8]. In addition, a clear heat transfer law within the fractures is of great significance to the establishment of EGS heat recovery models [9].

**Citation:** Gao, X.; Zhang, Y.; Hu, Z.; Huang, Y. Numerical Investigation on the Influence of Surface Flow Direction on the Heat Transfer Characteristics in a Granite Single Fracture. *Appl. Sci.* **2021**, *11*, 751. https://doi.org/10.3390/ app11020751


Received: 2 December 2020 Accepted: 11 January 2021 Published: 14 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Numerous studies have attempted to quantify the roughness of fracture surfaces. For example, Patton [10] proposed the concept of the wave angle and established the relationship between wave angle and fracture morphology. Mandelbrot [11] proposed the concept of fractal geometry and developed a fractal dimension framework to describe the geometrical characteristics of fracture rough surfaces. Barton [12] used the joint roughness coefficient (*JRC*) to characterize the section geometry of fracture surfaces to determine 10 fracture types with an approximate length of 10 cm. Xie [13] demonstrated the fractal dimension to be the measured values of the joint roughness coefficients.

Experimental and numerical simulations have demonstrated the fluid flow and heat transfer processes in fractures to be affected by many factors, including the aperture, roughness, type of fluid, injection flow rate, initial temperature of rock, and rock thermophysical parameters [14]. Early scholars used the plate model to describe the heat transfer characteristics in fracture. However, Tsang and Brown [15,16] pointed out that the parallel plate theory may cause flow estimation errors of 1–2 orders of magnitude. In an experimental study on the seepage heat transfer in a single fracture of granite, Zhao and Tso [17] revealed the single fracture with a certain roughness to have a better heat transfer performance than the smooth-plate fracture. Similarly, the two-dimensional coupled heat transfer model of a single fracture established by He [18] also demonstrated a strong correlation between the local convection heat transfer coefficient and the fracture profile. Neuville [19] established a self-affine rough fracture coupling model that assumes a reduced heat transfer efficiency in a rough fracture with constant temperature due to the channel effect. Research carried out by Huang [20] also showed that a larger relative roughness increases the flow friction and significantly reduces the heat exchange. Luo [21] investigated the influence of initial rock temperature on the heat transfer efficiency of fractures, demonstrating a positive relationship between the heat transfer rate and rock temperature. The two-dimensional single-fracture heat transfer model established by He [22] shows that the effect of fracture surface roughness on heat transfer intensity decreases with an increasing injection flow rate. Confining pressure has a significant impact on the fracture aperture. Shu [23] employed a novel experimental device to simulate the evolution of the hydraulic and heat transfer characteristics of fractures under a confining pressure. Results show that both the heat transfer rate and total heat transfer coefficient decrease with an increasing confining pressure. Zhang [24] established a three-dimensional numerical model based on threedimensional laser scanning to evaluate the effects of rock temperature, water flow velocity, roughness, and fracture aperture on the heat transfer coefficient. Simulation tests reveal the dominant influence of the flow rate on the rock roughness, followed by the aperture size. Andrade et al. [25] simulated the flow and heat transfer between two-dimensional parallel rough surfaces, demonstrating the minimal effect of wall roughness on the heat transfer for the characteristic parameter *Pe* = *u*(*d*/2)2/*L* < 200. Zhang [26] performed convection heat transfer tests of carbon dioxide in a single fracture to establish a numerical model and proposed the rough fracture channel effect. Experimental and numerical simulations can increase our understanding of the geothermal field model. Fox et al. [27] established a discrete fracture network and investigated the relationship between hydraulic aperture and fluid runoff, concluding that variations in the fracture aperture can cause flow channeling and reduce the heat transfer area. Furthermore, several studies [28,29] have employed Monte-Carlo [30,31] stochastic simulations to build fracture network models with different geometric parameters, and in particular, the EGS site model to predict the heat recovery performance. Bruel [32] integrated experimental data into a single model used to generate an estimation model of the main hydraulic parameters. Chen [33] established a geothermal reservoir model with a rough fracture surface based on small-scale research results, revealing the constant heat transfer coefficient (HTC) recommended in previous studies to underestimate the final outlet fluid temperature in the case of rough fractures. These studies demonstrate that research on the influence of fracture morphology on the fluid flow and heat transfer process at the experimental scale can be applied to the actual site model [28,34].

In summary, previous studies have shown that the geometry of the fracture surface has a profound influence on the fluid flow and heat transfer process in the fracture [35,36]. According to the different influencing factors, many scholars [23,37,38] have designed experiments and carried out corresponding numerical simulations to determine how these influence factors affect the heat transfer process in fractures. The results of such studies are very important for the development of EGS stimulation technology.

The aforementioned research generally focuses on the influence of fracture morphology on fluid heat transfer, while studies on the effect of the angle between the direction of the fluid flow and the morphology of the fracture wall on the heat transfer performance are limited. In addition, because of the geological history, the fracture aperture of EGS reservoirs range from the micron to centimeter scale [39]. This makes it difficult to obtain the temperature distribution of rock fractures and thermal mediums under laboratory conditions. Hence, in the current study, we established a single fracture heat transfer model with a random geometry profile that was subsequently verified by experiments. Four cases with fracture profiles and varying angles between flow direction were set up to simulate and explore the heat transfer performance of distilled water through fractures. The influence of the fluid flow direction on the fracture heat transfer characteristics was then discussed.

#### **2. Numerical Simulations**

#### *2.1. Mathematical and Physical Models*

A three-dimensional single fracture model was established to study the influence of flow direction on the heat transfer characteristics in a granite single fracture. The geometry profile of the fracture was constructed using a random method (Figure 1a), while the rough fracture wall was built by stretching in the vertical direction (Figure 1b,c).

**Figure 1.** Fracture surface and surface flow direction distribution. (**a**) 2D random fracture profile; (**b**) pseudo-threedimensional fractured rough surface stretched towards the page direction. (**c**) 3D view of pseudo-three-dimensional fracture rough surface.

The surface flow of the fracture surface was assumed to have four different directions, with the 0◦ direction taken as the *z*-axis direction (Figure 1b). The angle between the surface flow direction and the positive direction of the *z*-axis is denoted as α. The surface flow direction had a gradient of 30◦ and was rotated counterclockwise around the *z*-axis by 30◦, 60◦, and 90◦ to obtain the fracture wall. The rough single-fracture numerical model shown in Figure 2 was obtained by embedding these fracture walls into a cylindrical model with a size of 100 × 50 mm (length × diameter). The vertical displacement of the fracture in the model was determined as 0.3 mm by translation method. Table 1 lists the geometric parameters of the four fracture models. The surface area and inlet area

deviations of the four fractures were 1.87% and 0.33%, respectively. Therefore, at the scale of this study, the effects on the simulation results due to differences in surface area and inlet area are negligible.

**Figure 2.** Fracture seepage heat transfer model.


**Table 1.** Geometric parameters of the fracture.

For natural fracture, the Reynolds number is much smaller than the critical Reynolds number (Re < 1800 [35]) due to its small aperture and low flow rate of underground fluid. Thus, the fluid flow in a granite fracture remains laminar. Convective heat transfer in three-dimensional rock fractures was modeled using FLUENT 6.3 [40] by employing a finite volume method (FVM) to solve the Navier–Stokes equations.

The governing equations for the conservation of momentum, mass, and energy in the fracture are described as follows:

Continuity equation:

$$\frac{\partial \rho\_w}{\partial t} + \nabla \cdot (\rho\_w \mathbf{V}) = 0 \tag{1}$$

Momentum equation:

$$\frac{\partial(\rho\_w V)}{\partial t} + \nabla(\rho\_w V V) = -\nabla P + \nabla \cdot \left(\mu \left[\left(\nabla V + \nabla V^T\right) - \frac{2}{3}\nabla \cdot V I\right]\right) \tag{2}$$

Energy equation:

$$\frac{\partial}{\partial t}(\rho\_w E) + \nabla \cdot (\mathcal{V}(\rho\_w E + p)) = \nabla \cdot \left(\lambda\_w \nabla T + \mu \left[\left(\nabla V + \nabla V^T\right) - \frac{2}{3} \nabla \cdot \mathbf{V} I\right] \cdot \mathbf{V}\right) \tag{3}$$

where *ρ<sup>w</sup>* (kg/m3) and *ρ<sup>r</sup>* (kg/m3) are the densities of water and rock, respectively; *V* (m/s) is the velocity vector; *p* (Pa) is the pressure; *I* is the unit tensor; *E* is the total energy; *λ<sup>w</sup>* (W/(m· ◦C)) is the fluid thermal conductivity; *μ* (Pa s) is the fluid viscosity.

In EGS frameworks, water is typically used as a heat exchange medium. However, the physical properties of water are more susceptible to temperature than pressure. Hence, the relationship between the physical properties of water and temperature was described by the following empirical formula [41]:

Dynamic viscosity of water:

*<sup>μ</sup><sup>w</sup>* = 1.3799 − 0.0212*<sup>T</sup>* + 1.3604 × <sup>10</sup>−4*T*<sup>2</sup> − 4.6454 × <sup>10</sup>−7*T*<sup>3</sup> + 8.9043 × <sup>10</sup>−10*T*<sup>4</sup> − 9.0791 × <sup>10</sup>−13*T*<sup>5</sup> <sup>+</sup>3.8457 <sup>×</sup> <sup>10</sup>−16*T*<sup>6</sup> (4)

Specific heat capacity of water:

$$C\_{p,w} = 12010 - 80.4 \times T + 0.3 \times T^2 - 5.4 \times 10^{-4} T^3 + 3.6 \times 10^{-7} T^4 \tag{5}$$

Thermal conductivity of water:

$$
\lambda\_w = -0.8691 + 8.9 \times 10^{-3} T - 1.5837 \times 10^{-5} T^2 + 7.9754 \times 10^{-9} T^3 \tag{6}
$$

Water density:

$$
\rho\_w = 838.4661 + 1.4005 \times T - 3 \times 10^{-3} T^2 + 3.7182 \times 10^{-7} T^3 \tag{7}
$$

where *T* ( ◦C) is the water temperature and which ranges from 10 to 100 ◦C.

The fracture inlet and outlet employ mass flow rate and outflow boundary conditions, respectively. The outer surface of the rock is considered to be an adiabatic boundary, while the inner wall of the fracture is considered as an impervious and no-slip boundary.

#### *2.2. Model Meshing and Mesh Independence Verification*

In order to focus on fluid flow and heat transfer in fractures, we employed structured hexahedral grids in fracture spaces (Figure 3a) and an increased grid density around fractures (Figure 3b,c), while an unstructured tetrahedral grid was used for the computational domain outside the fracture. Furthermore, we selected α = 0◦ as the verification sample to verify that the simulation results are not affected by the number of meshes. Table 2 lists the four grid numbers for α = 0◦ and Figure 3 presents the grid characteristics of the fourth grid number.

**Figure 3.** Mesh features of mesh 4. (**a**) Fracture surface mesh characteristics; (**b**) overall mesh distribution; (**c**) mesh distribution characteristics around the fracture.

**Table 2.** Four grid numbers for α = 0◦.


Under the same simulation conditions, four models with different numbers of grids were constructed and the average temperature of the fluid at the outlet was monitored (Figure 4). Mesh case 1 is distorted due to the presence of grids with a high value and aspect ratio. When the number of meshes exceeds 1.3 million, the outlet temperature of different mesh cases is generally consistent with time. In order to accurately reflect the geometry of the fracture surface and save computer resources, we selected mesh case 3 for further numerical calculations. Table 3 details the grid numbers for the four models.

**Figure 4.** Variation of outlet temperature with time for the four mesh cases.


**Table 3.** Number of simulated sample grids for the four models.

#### *2.3. Initial Conditions and Simulations*

The initial temperature of the rock was set to 90 and 70 ◦C to prevent a water phase transformation. The initial flow rates were 10, 20, 30, and 40 mL/min, respectively, depending on the performance of the plunger pump. The injection temperature of water was 25 ◦C without considering the change of room temperature. Table 4 lists the simulation conditions.

**Table 4.** Numerical simulation conditions.


Calculations were performed following the unsteady double precision method and the second-order implicit mode to improve calculation accuracy. The SIMPLEC (semiimplicit method for pressure-linked equations consistent) algorithm, an improved SIMPLE algorithm, was adopted to amend the fluid pressure and flow rate, while the PRESTO (PREssure STaggering Option) scheme was used for the pressure discretization. The second order discretization scheme was used for the convection terms. The residual convergence criteria for the computational equations was set to 10−<sup>5</sup> and the time step size for the transient simulations time was equal to 1 s. The total solution time was taken as 30 min.

#### **3. Model Verification**

Although the establishment of a realistic fracture model proves to be a complicated task, a single-fracture granite model with a flat and smooth surface can validate the proposed model. Thus, experimental and numerical simulation results for smooth-plate fractures are compared.

A fracture heat transfer laboratory simulation system was independently developed by our research group for the experiments. Figure 5 depicts the experimental system, which can be divided into five key components including the seepage, confining pressure, heating, core holder, and data measurement systems. Six temperature sensor mounting holes were preset on the top of the holder to ensure direct contact between the sensor and the outer surface of the rock sample. The data measurement system monitored the rock and fluid temperatures on the outer surface of the core and at the core outlet in real time, respectively.

**Figure 5.** The fracture heat transfer laboratory simulation system. (**a**) Core holder system; (**b**) a granite sample with a smooth-plate single fracture in the holder.

The rock sample with smooth-plate fractures shown in Figure 6a,b was tested to verify the correctness of the numerical model, while a numerical model (Figure 6c) was established to simulate the test results. The rocks used in this study were obtained from important geothermal resource targets in the northern Songliao Basin of China. The granite was processed into a cylindrical specimen with a diameter and height of 50 and 100 mm, respectively. The thermophysical properties of the granite were obtained from tests. The specific heat capacity measuring device (BBR series, made in Xiangtan Xiangyi Instrument Co., Ltd., Xiangtan, China) based on cooling method and the thermal conductivity scanner (TCS, made in Lippmann and Rauen GbR, Schaufling, Germany) based on optical scanning principle are used to measure the specific heat capacity and thermal conductivity of rock, respectively. Figure 7 presents the variations in specific heat capacity and the thermal conductivity of rock with temperature.

**Figure 6.** Smooth-plate fracture and mesh model. (**a**) Temperature sensor arrangement. (**b**) Flow direction. (**c**) Mesh generation characteristics of numerical model.

**Figure 7.** Variations in thermal conductivity and specific heat capacity of rock with temperature.

The experimental flow rates were set as 5, 10, 20, and 30 mL/min for temperatures of 60, 70, 80, and 90 ◦C, respectively. Both the boundary and initial conditions followed those detailed in Section 2.3. Numerical simulations were performed based on the experiments to ensure that the results were not affected by other factors. The unsteady numerical simulation of fluid temperature *Toutlet* at the fracture outlet for a given time was calculated as follows:

$$T\_{outlet} = \frac{\int\_{A\_{\rm out}} T\_w(x, y, z) \mu\_w(x, y, z) \rho\_w(x, y, z) c\_{p, w}(x, y, z) dA}{\int\_{A\_{\rm out}} \mu\_w(x, y, z) \rho\_w(x, y, z) c\_{p, w}(x, y, z) dA} \tag{8}$$

where *Aout* (m2) is the area of the fracture outlet.

Figure 8 clearly describes the comparison between the numerical simulation and experimental results for a smooth single fracture. Although there are some errors in the outlet temperature and the rock surface temperature obtained from the experiments and numerical simulation, the maximum errors are only 0.8% and 1.3%, respectively. This indicates that the numerical model established using FLUENT 6.3 can accurately reflect the convection heat of single-fracture granite.

**Figure 8.** Comparison of experimental and numerical results (*t* = 10 min). (**Left**) Outlet temperature varies with injection flow rate at different initial rock temperatures. (**Right**) Surface temperature varies with sensor position at different injection velocities.

#### **4. Data Processing**

#### *4.1. Heat Extraction Rate*

The heat extraction rate can measure the heat absorbed by water flowing through a rock fracture per unit time and is calculated as follows [26]:

$$\mathcal{Q} = (H\_{\text{outlet}} - H\_{\text{inlet}})M\_\prime \tag{9}$$

where *M* (kg/s) is the flow rate; *Hinlet* (kJ/kg) and *H*outlet (kJ/kg) are the water enthalpy at the inlet and outlet, respectively. The enthalpy of water can be determined by numerical simulations of the inlet and outlet temperature and pressure.

#### *4.2. Total Heat Transfer Coefficient*

For fracture channels with complex surface morphology, the total heat transfer coefficient can describe the overall convective heat transfer performance within the fracture. The total heat transfer coefficient can be calculated from Newton's law of cooling [42]:

$$h\_{total} = \frac{q\_f}{T\_f - T\_w} \,\tag{10}$$

where *q <sup>f</sup>* (W/m2) is the average heat flux on the inner wall of the fracture; *Tf* ( ◦C) is the average temperature on the inner wall of the fracture; *Tw* ( ◦C) is the average temperature of the water in the fracture. These parameters can be calculated from the numerical simulation results as follows:

$$q\_f = \frac{\int\_{A\_f} q(\mathbf{x}, y, z) \mathbf{d}A}{\int\_{A\_f} \mathbf{d}A},\tag{11}$$

$$T\_f = \frac{\int\_{A\_f} T(x, y, z) \mathrm{d}A}{\int\_{A\_f} \mathrm{d}A},\tag{12}$$

$$T\_{\rm{IV}} = \frac{\int\_{V\_f} T(\mathbf{x}, \mathbf{y}, z) \rho(\mathbf{x}, \mathbf{y}, z) \mathbb{C}\_p(\mathbf{x}, \mathbf{y}, z) u(\mathbf{x}, \mathbf{y}, z) \mathrm{d}V}{\int\_{V\_f} \rho(\mathbf{x}, \mathbf{y}, z) \mathbb{C}\_p(\mathbf{x}, \mathbf{y}, z) u(\mathbf{x}, \mathbf{y}, z) \mathrm{d}V}, \tag{13}$$

where *Af* (m2) is the surface area of the fracture; V*<sup>f</sup>* (m3) is the volume of the fracture.

#### *4.3. Local Heat Transfer Coefficient*

The local heat transfer coefficient can describe the local heat transfer characteristics in fractures and is an important parameter reflecting the heat transfer characteristics in the flow direction. Here, we calculated the discrete form of the local heat transfer coefficient via Equation (10) as follows:

$$h\_{\mathbf{x}} = \frac{q\_{f\mathbf{x}}}{T\_{f\mathbf{x}} - T\_{\mathbf{xx}}\prime} \tag{14}$$

where *q f x* (W/m2), *Tf x* ( ◦C), and *Twx* ( ◦C) are the heat flux, temperature, and fluid temperature at position *x* of the fracture wall, respectively. Furthermore, the Reynolds number (*Re*) can be calculated using the following formula for seepage in single fracture of rock:

$$R\_{\mathfrak{k}} = \frac{2M}{d\mu\_w},\tag{15}$$

where *M* (kg/s) is the mass flow rate and *d* (m) is the fracture aperture.

The overall Nusselt number (*Nu*) can be determined as follows:

$$N\_{\mu} = \frac{h\_{total} d\_e}{\lambda\_w} \tag{16}$$

where *de* (m) is the equivalent hydraulic diameter of a channel with arbitrary shape.

#### **5. Results and Discussion**

#### *5.1. Temperature Field Distribution*

Table 5 lists the outlet water temperature determined from the 32 simulations. The outlet water temperature is maximized at α = 90◦, followed by models α = 0◦, α = 30◦, and α = 60◦.

**Table 5.** Results from 10-min numerical simulation.


Figure 9 depicts the distribution characteristics of the temperature fields on the fracture surfaces of four different model cases at several injection flow rates. With the injection of low temperature water, the temperature inside the fracture begins to decrease and the cold front pushes towards the outlet. This phenomenon is consistent with previous studies [37,43,44]. The temperature at the center of the fracture typically exceeds that on both sides of the fracture. This is attributed to the greater thermal replenishment at the center of the fracture by the cylinder model compared to the two sides. Increasing the flow rate enhances the rate of the temperature reduction in the fracture. The temperature on the fractured surface of model α = 0◦ is distributed symmetrically along the *z*-axis due to the symmetry of the model. At α = 30◦ and α = 60◦, the fracture surface temperature in the positive direction of the *x*-axis exceeds that in the negative direction. This indicates that the angle between the fracture morphology and flow direction has an important influence on the convection heat transfer process.

Figure 10 presents the local temperature distribution along the *z*-axis at *x* = 0.02 m and *x* = −0.02 m. At α = 30◦ and α = 60◦, the local wall temperature at *x* = 0.02 m is lower than that at *x* = −0.02 m. This is related to the migration of the channel formed by the fracture profile towards the *x*-axis, which makes the fracture geometry profile no longer perpendicular to the flow direction. The fluid preferentially flows through the fractures in an easily accessible path where the flow direction is deflected. In addition, Figure 10 shows that the temperature curve at low flow rates (10 mL/min) reaches its peak value, while the temperature curve at high flow rates (30 mL/min) does not. This indicates that the cryogenic fluid injected at high flow rates penetrates the fractures more quickly and results in a thermal breakthrough in a short time, which causes rapid wall temperature reduction and no timely thermal replenishment. Fluids with low flow rates are more likely to utilize matrix heat transfer to compensate for the decrease in wall temperature caused by convection heat transfer.

Figure 11 illustrates the streamline distribution characteristics of the models with α = 30◦ and α = 60◦. The flow direction is observed to slightly deflect in the fracture towards the *x*-axis, which results in the temperature difference shown in Figures 9 and 10. This deflection is attributed to the presence of a gently inclined microfractured surface in front of the flow field, preventing the flow from passing. Figure 12 simplifies fracture surface morphology to account for the deflection of the fluid flow direction. A fractured surface consists of a limited number of microfractured surfaces. When the direction of the fluid injection is not parallel to the microfracture dip, the fluid will slightly deflect towards a blunt angle between the flow direction and the microfracture surface strike.

**Figure 9.** Temperature field distribution on the fracture surface (*t* = 10 min).

**Figure 10.** Local temperature distribution along the flow direction at different locations of the *x*-axis on the inner wall of fracture (solid line: *x* = −0.02 m, dashed line: *x* = 0.02 m, *t* = 2 min).

**Figure 11.** Streamline distribution characteristics of models with α = 30◦ and α = 60◦ (*V*<sup>0</sup> = 10 mL/min).

**Figure 12.** Schematic illustration of fluid flow direction deviation. (**a**) The up- and down-hill direction of the fluid migration. (**b**) Top view: solid and dotted lines are the ideal flow direction path and the flow deflection path, respectively.

#### *5.2. Heat Extraction Rate*

The average Reynolds and Nusselt numbers in the fracture are calculated based on Equations (15) and (16), respectively (Figure 13). Based on the simulation conditions, the Reynolds number calculated using Equation (16) ranges from 14.5 to 53.3, which is much lower than the critical Reynolds number leading to the change of laminar flow state [35]. The results of Reynolds number calculation also show that the assumption of using laminar flow model in Section 2.1 is reasonable. The difference of thermophysical properties of water results in a greater Reynolds number at initial temperature *T*<sup>0</sup> = 90 ◦C compared to that with initial temperature *T*<sup>0</sup> = 70 ◦C. The maximum Nusselt number at α = 90◦ is 15.8 (Figure 13b).

**Figure 13.** Reynolds number in relation to (**a**) the flow rate and (**b**) Nusselt number.

Considering the small scale of the model in this study, we used the cumulative heat extraction rate to represent the heat absorbed by water from the model. The cumulative heat extraction rate can be calculated as follows:

$$Q\_{\text{sum}} = \sum (H\_{\text{outlet},t} - H\_{\text{inlet},t})M.\tag{17}$$

Figure 14 presents the cumulative heat extraction rate *Qsum* (W) under different operating conditions. The cumulative heat extraction rate is maximized at α = 90◦, followed by α = 0◦, α = 30◦, and α = 60◦. Within the local range, increasing the roughness of the fracture does not necessarily enhance its heat transfer capacity. This is inconsistent with the results of other studies because other researchers use two-dimensional models [22,45]. Moreover, although the cumulative heat extraction increases with the injection flow rate under the same model [24,43], the increasing rate (Δ*Qsum*, see Figure 14a) will gradually decrease. Thus, blindly increasing the injection flow rate during EGS operation may be detrimental to heat recovery. At the same flow rate level, an increase of the initial temperature from 70 to 90 ◦C will enhance the cumulative heat extraction rate by approximately 45%. Figure 15 demonstrates the variation of heat extraction rate with time, indicating that a larger initial flow rate at the initial stage of the heat transfer will improve the heat extraction capacity. The heat loss on the fracture wall caused by high convection heat exchange caused by high flow rate is much greater than that supplied by the matrix by heat conduction, which causes a sharp decrease in heat extraction rate. The higher the flow rate, the faster the descent rate.

**Figure 14.** Histogram of cumulative heat extraction rate determined for 30 min.

**Figure 15.** Variation of heat extraction rate with time at several initial flow rates and initial temperatures.

#### *5.3. Total Heat Transfer Coefficient*

Figure 16 shows the variation of the total convective heat transfer coefficient with time. The total heat transfer coefficient curve can be divided into increasing and decreasing stages. During the increasing stage, the low temperature water injection causes severe convection heat exchange on the fracture wall, which can prevent cold water from penetrating the whole fracture. Therefore, the total heat transfer coefficient increases rapidly with the injection of fluid until it reaches its peak value. At the decreasing stage, the low temperature water enters the fractures continuously and the heating effect of the rock matrix begins to weaken with the decrease of temperature. In addition, the fluid and fracture wall temperatures further decrease, and the low temperature water begins to penetrate the fractures. Although increasing the initial flow rate and temperature can amplify the total heat transfer coefficient, increasing the flow rate will shorten the time required by the total heat transfer coefficient to reach its peak value. This is obviously not conducive to the sustainability of heat recovery. Figure 16 shows that the total heat transfer coefficient in descending order is α = 90◦, α = 0◦, α = 30◦, and α = 60◦.

**Figure 16.** Variations in the total heat transfer coefficient with time for various initial conditions. Top: *T*<sup>0</sup> = 90 ◦C. Below: *V*<sup>0</sup> = 10 mL/min.

#### *5.4. Local Heat Transfer Coefficient*

Figure 17 depicts the local heat transfer coefficients of the rock fracture surfaces. The local convection heat transfer coefficient is observed to gradually decrease along the outlet direction, which is consistent with the previous study [22,42]. This is attributed to the gradual decrease in the temperature difference between the fluid and fracture wall along the outlet direction. The fracture model shows that the fracture roughness of α = 30◦ and α = 60◦ are different in the three axes, while that of α = 0◦ and α = 90◦ are different in the two axes. The flow direction on the fractured surface is slightly deflection, and consequently the local convection heat transfer coefficient in the negative direction of the *x*-axis exceeds that of the positive direction (Figures 11 and 12). In addition, the distribution of the local heat transfer coefficient is strongly correlated with the rough fracture profile [45]. Figure 18 shows the distribution of the local heat transfer coefficients on the fracture profile at *x* = 0 m. Figure 18 shows that the model of fracture roughness from large to small at *x* = 0 m is α = 0◦, α = 30◦, α = 60◦, and α = 90◦. The fluctuation of local convection heat transfer coefficient increases with the increase of fracture roughness. At α = 90◦, the local heat transfer coefficients exhibit a logarithmic distribution in the smooth section along the *z*-axis. Similar to the fracture profile, the distribution curve of the local heat transfer

coefficient is characterized by peaks and troughs. The local heat transfer coefficient curves at α = 0◦, α = 30◦, and α = 60◦ reveal that the local heat transfer coefficient at the peaks of the fracture profile is consistently greater than that at the troughs. Moreover, increasing the injection flow rate not only improves the local heat transfer coefficient, but also has an enhancing effect that makes the difference between the peaks and troughs of the local heat transfer coefficient more evident. The enhanced effect can be observed via the increased slope of the local heat transfer coefficient curve. Figure 19 shows the enhancement effect of increasing the flow rate on the convective heat transfer intensity.

**Figure 17.** Distribution of local heat transfer coefficients on rock fractures as water flows through the fracture (*V*<sup>0</sup> = 10 mL/min. *T*<sup>0</sup> = 90 ◦C. *t* = 10 min).

**Figure 18.** *Cont.*

**Figure 18.** Distribution of local heat transfer coefficients along the outlet direction at different flow rates.

**Figure 19.** Variation of *hx* with the flow rate at adjacent peak and troughs of the fracture profile.

Figure 20 presents the distribution of the local heat transfer coefficient, wall temperature, and average flow rate on the fracture profile at *z* = 0.05 m. The curves of the average flow rate and local heat exchange coefficient become increasingly complex as α increases from 0◦ to 90◦. Since the fracture model is constructed by translating along the *y*-axis, the distribution of fracture aperture in each direction of the model is not uniform. Therefore, a position with a higher absolute gradient results in a smaller fracture aperture, which results in a low flow rate as the water flows through the path of the small aperture. In contrast, fractures with larger aperture will form a dominant channel where the flow rate is greater. Compared with the temperature, the difference of flow rates at fractures with different apertures is the main factor causing uneven distribution of local heat transfer coefficient. In addition, a strong positive correlation is observed between the local heat exchange coefficient and the flow rate distribution, which indicates that the flow rate is the main factor affecting the heat transfer intensity.

**Figure 20.** Local heat transfer coefficient, wall temperature, and average velocity distributed along the *x*-direction at *z* = 0.05 m (*V*<sup>0</sup> = 20 mL/min, *T*<sup>0</sup> = 90 ◦C).

#### *5.5. Limitations of the Present Study*

The rough fracture surface used in this study is actually a pseudo-3D fracture surface that ignores the anisotropy of roughness to a certain extent. Therefore, models for this study are not common in actual natural environments. Moreover, the scale of the fracture model may weaken the impact of the proposed flow direction deflection on the heat transfer effect. The heat transfer characteristics of inhomogeneous three-dimensional rough fractures need to be further studied. However, the focus of this work is to study how different angles between flow direction and fracture profile affect the heat transfer process in fractures. In fact, a certain size fracture surface is composed of a limited number of microfractures with occurrence, which have different angles with the flow direction of fluid. From this point of view, this work can actually provide basic knowledge for EGS research.

#### **6. Conclusions**

In this study, a hydro-thermal coupling model was developed to investigate the influence of fluid flow direction on the heat transfer characteristics of fractures. The accuracy of the proposed numerical model was verified by laminar convection heat transfer tests of smooth-plate fractured rock samples. A detailed analysis of the temperature field distribution, total heat transfer coefficient, and local heat transfer coefficient of the model were provided based on the numerical simulations. The main conclusions are as follows:

(1) At α = 30◦ and α = 60◦, the temperature of the fracture surface in the negative direction of the *x*-axis is higher than that in the positive direction. This is because the flow direction of the fluid is not parallel to the orientation of the microfracture surface that prevents the fluid from flowing forward, resulting in fluid deflection.

(2) When the initial temperature is raised from 70 to 90 ◦C, the cumulative heat extraction rate *Qsum* increases the yield by 45%. Although the cumulative heat extraction increases with the injection flow rate, the rate of increase gradually decreases. An increase in the injection flow rate from 10 to 40 mL/min enhances the peak heat extraction rate. However, the increase in flow rate speeds up the rate of reduction of the heat extraction rate.

(3) Variations in the total heat transfer coefficient with time can be divided into increasing and decreasing stages. The total heat transfer coefficients decrease with α. This indicates that the presence of rough fractured surfaces in the forward direction of the fluid does not necessarily increase the convection heat transfer intensity between the fluid and the fractured surfaces. Increasing the flow rate and temperature can enhance the total heat transfer coefficient. However, a higher flow rate will shorten the time required by the total heat transfer coefficient to reach its peak value, while also increasing the rate of reduction of the total heat transfer coefficient in the decreasing stage.

(4) The local heat transfer coefficient gradually decreases along the direction of the water flow, which is consistent with the conclusions drawn by Bai [46]. An increasing flow rate not only increases the local heat transfer coefficient, but also has an enhancing effect that highlights the influence of the fracture profile on the local heat transfer coefficient distribution.

(5) There is a consistently strong correlation between the local heat transfer coefficient and the fracture profile, independent of the axis. In particular, although a change in the angle between the flow direction and the fracture profile can affect the local heat transfer coefficient, it does not affect the distribution characteristics of the local heat transfer coefficient.

**Author Contributions:** Conceptualization, writing—original draft preparation, X.G.; supervision and review, Y.Z.; formal analysis, data curation, Z.H.; software, Y.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Natural Science Foundation of China (Nos. 41772238), by the New Energy Program of Jilin Province (Nos. SXGJSF2017-5), and the Interdisciplinary Fund of Jilin University (Nos. 101832020DJX072).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used to support the findings of this study are available from the first author or the corresponding author upon request.

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

#### **References**


## *Article* **Well-Doublets: A First-Order Assessment of Geothermal SedHeat Systems**

**SeyedBijan Mahbaz 1,2, Ali Yaghoubi 1, Alireza Dehghani-Sanij 2,\*, Erfan Sarvaramini 3, Yuri Leonenko 1,2,4 and Maurice B. Dusseault 1,2**

	- Waterloo, ON N2L 3G1, Canada

**Abstract:** Renewable and sustainable energy sources can play an important role in meeting the world's energy needs and also in addressing environmental challenges such as global warming and climate change. Geothermal well-doublet systems can produce both electrical and thermal energy through extracting heat from hot-water aquifers. In this paper, we examine some potential challenges associated with the operation of well-doublet systems, including heat conductivity, chemical, and mechanical issues. In these systems, geomechanics issues such as thermal short-circuiting and induced seismicity arise from temperature and pressure change impacts on the stress state in stiff rocks and fluid flow in fractured rock masses. Coupled chemical processes also can cause fluid channeling or formation and tubular goods plugging (scaling) with precipitates. Mechanical and chemical disequilibrium conditions lead to increased production uncertainties; hence risk, and therefore coupled geo-risk assessments and optimization analyses are needed for comparative commercialization evaluations among different sites. The challenges related to heat transfer processes are also examined. These studies can help better understand the issues that may arise during the operation of geothermal well-doublet systems and improve their effectiveness, subsequently reducing associated costs and risks.

**Keywords:** geothermal; well-doublet system; sustainability; disequilibrium; thermomechanical effects; chemical coupling; climate change

#### **1. Introduction**

Energy is an essential and unavoidable need in today's world and has played a crucial role in human civilization [1–3]; demand is increasing, with the principal reasons being the world's growing population and the general desire for an improved life quality [4–7]. Fossil fuels comprise over 80% of primary energy sources [8,9], but effects arising from greenhouse gas (GHG) emissions and air pollution, leading to environmental degradation, global warming, and climate change, are of increasing concern [10–13]. Geothermal energy, a sustainable and green source, is now meeting some communities' thermal and electrical energy requirements [10,11,14]; it can be reliable in the long term, is environmentally friendly [7,15], and could help reduce fossil fuel consumption [8,9].

#### *1.1. Geothermal Energy Characteristics*

Geothermal energy can be categorized as volcanic, sedimentary, geo-pressured, hot dry rock (HDR) or enhanced (engineered) geothermal (EG), or shallow [7,16]. It should be noted that shallow geothermal energy can be defined in terms of shallow geothermal systems

**Citation:** Mahbaz, S.B.; Yaghoubi, A.; Dehghani-Sanij, A.R.; Sarvaramini, E.; Leonenko, Y.; Dusseault, M.B. Well-Doublets: A First-Order Assessment of Geothermal SedHeat Systems. *Appl. Sci.* **2021**, *11*, 697. https://doi.org/10.3390/app11020697

Received: 28 December 2020 Accepted: 8 January 2021 Published: 13 January 2021


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

based on ground-source heat pumps that are installed at depths less than 500 m [17]. Geothermal resources can also be divided into three groups based on reservoir fluid temperatures [18]: low, medium, and high enthalpy resources. Their common point is accessing natural underground heat sources, whether from dry steam reservoirs (electrical power) to the normal warmth of deeply buried rocks and their interstitial fluids (district heating). Although such energy is often classified as renewable because of the continuous heat flux in the Earth's crust, realistically, the heat renewal rate (≈100–1000 years) is far slower than the design life-span of a commercial project (25–40 years) unless strong hot fluid recharge is occurring. The extraction system and energy rate must be compatible with project life needs [19], and in most cases, from an engineering project life-span definition, geothermal energy is not renewable. In other words, the life-span of any commercial geothermal project is strongly dependent on the geological characteristics of the resources and the appropriate sizing of the power/heat plant.

Once a spatially limited volume of the reservoir has been "heat mined," the heat recharge rate from crustal thermal conduction is too slow to sustain commercial use unless an external heat source is being used to "recharge" the reservoir heat content. However, heat is everywhere available at depth, so the project design basis should include continued access to previously unmined heat (i.e., drilling more wells on an ongoing basis or recharging the repository heat).

Economically, the amount of energy available and the rate of energy production are essential criteria for project performance. Although initial reservoir temperature (To) is the dominant variable, which defines the energy content, the amount of energy production is also related to the reservoir pressure—P (natural or maintained if fluid injection is being exploited). Overall, a higher temperature and pressure reservoir is better [8], but there are a variety of other variables such as reservoir volume, permeability, porosity, storativity, and natural recharge rate (if any) that are important to the commercially viable project duration. Site assessment factors include reservoir geometry (depth, thickness, and lateral extent), flow characteristics (permeability heterogeneity, fracture vs. porous media flow, channels, etc.), the potential for natural recharge, the tectonic state and stress field, and structural geology factors [20].

Whether a reservoir is open or closed to natural high-temperature fluid influx impacts life-span assessment: in a closed system with a constant production and no injection, reservoir pressure decreases continuously, requiring aggressive pumping or leading to a rapid loss of water production, although the rock mass may remain at high temperature. These systems may be considered for re-injection operations, as a closed sedimentary geothermal formation can be developed with a well-doublet (or multi-well) system. In contrast, an open system has a more stable pressure but may require flow rate management and assessment of the recharge rate and energy content (production temperature Tin) to assure commercial viability [19]. An open system with an "infinite" flow capacity can be exploited with single pumped production wells, but the geothermal fluids have to be disposed of in the subsurface into suitable (likely shallower) saline aquifers.

In the current paper, we use the increasingly common term "SedHeat" (sedimentary heat) for systems where warm and hot fluids are being exploited from sedimentary rocks of sufficient permeability and porosity at depth [21,22], in contrast to steam systems (wet or dry), or to geothermal heat found in low-permeability rock masses (granites, shales, tight volcanic or metamorphic rocks, etc. [23,24]).

#### *1.2. The Geothermal Doublet*

The well-doublet design (Figure 1) is the simplest approach to energy extraction in a closed or limited recharge SedHeat system, whether it is a deep liquid reservoir or a shallow ground-source heat pump system using an unconfined aquifer. From one well equipped with a lifting pump, the SedHeat aquifer supplies hot-water from which energy is extracted by a heat exchanger to be sent for power generation or direct heat use. Cooled geofluid re-injection takes place into another well placed at an optimum distance designed to meet the project's desired power profile (power = electrical and heat energy rate) [25] and sustain flow. The maximum energy rate (work or power) available, ignoring all of the minor hydraulic losses (pumps, fitting, etc.) and irrecoverable heat losses, is

$$
\dot{\mathbf{E}} = \Delta \mathbf{T} \cdot \dot{\mathbf{Q}} \cdot \mathbf{c\_P} \tag{1}
$$

where . <sup>E</sup> is the energy rate (power) in J/s, <sup>Δ</sup>T = (Tin –Tout) in K, . *Q* is the flow rate in kg/s, and cp is the fluid heat capacity in J/kg·K. Tin and Tout (units of K) refer to the inlet and exit temperatures in the energy (electricity + heat) extraction system, likely a heat exchanger. The fluid involved may be the natural fluids in the case of SedHeat development, or in the case of EG systems, it is likely to be an introduced fluid such as water or supercritical CO2 [26].

**Figure 1.** Schematic illustration of a typical shallow geothermal well-doublet system with a single primary loop, no heat exchanger, in an unconfined aquifer (modified from [27]).

Specifically, Tin is not the original reservoir To; it is the temperature at the surface going into the energy system (the heat exchanger between the primary loop and the secondary power/heat system loop). It will eventually change with time as cool injected fluids begin to impinge upon the production well. Nor is Tout equal to the ambient temperature; Tout will always be higher because of thermodynamic constraints. Unextracted heat in the primary loop fluid is returned to the reservoir through fluid re-injection; if there are low-grade waste heat sources to increase Tout through a downstream heat exchanger, the project life may be extended, and the low-grade heat exploited. The parameter . *Q* is mainly a function of the aquifer's permeability, thickness, and layering, the amount of energy desired or achievable, the aquifer's thermal properties and heat recovery time [8], and several secondary factors. A reasonable value for cp is 4180 J/kg·K for pure freshwater, somewhat different for geofluids of different salinities. For fluids such as supercritical CO2, different specific heats, fluid densities, and viscosities must be used, and these may be quite sensitive to P and T, leading to additional non-linear effects in modeling heat extraction.

Using a recirculating well-doublet in a binary geothermal system means that remnant useable heat in the geofluid (e.g., 50–70 ◦C) is reinjected after electrical power generation, but if space heating is needed, additional heat can be stripped from the fluids, lowering re-injection T. Furthermore, coupling a deep doublet with a shallow ground-source heat pump system allows shallow seasonal heat storage and recovery (a thermal repository) and leads to higher efficiency and lower fossil-fuel use where there is a cooling deficiency in the summer [28], as in northern climates. Using a deep multi-well system allows injection/production strategies to be devised to maximize the system's value, while also exploiting the shallow heat repository in the winter when high heat demand exists. The

design of such a system depends on the estimated power/thermal energy profiles for the year, combined with a safety factor, and noting that design in Canada (for example) must be to meet November–March energy needs.

If natural or planned thermal recharge is rapid enough, considering the recovery time of the system and the injection/production strategy, a geothermal SedHeat system might approach being a renewable and sustainable energy source [25]. In other words, the deep system life-span can be extended by designing for slow energy extraction relative to the system volume, re-injection of unused heat, extracting value from sources of waste heat, and even using seasonal solar and excess wind power to recharge the thermal repository. It is possible to prolong the system's energy provision profile by keeping heat extraction at a sustainable level or by designing it initially so that the combined electrical power and heat outputs extend the long-term utility of the project. In addition, it may be possible for a SedHeat system to generate power and heat for 20–30 years, but in the later life stages when the Tin is degraded, and power generation is suffering, the original wells may be operated for heat production, but new sources developed for power generation. In this scenario, given that the capital investment has been recovered, the initial wellbore system can be economically continued for heat provision only. Geothermal energy developed in this way is inherently sustainable because of the huge masses of warm rock in the crust and within reasonable drilling depth. SedHeat is, of course, limited to sedimentary basins, but other technologies are being explored for low permeability rock masses such as crystalline rocks.

#### *1.3. Disequilibrium Processes*

Leaving aside initial geological conditions, the causes of disequilibrium that result in an inadequate or a drop in E are based on mechanical, geochemical, and heat conductivity ˙ conditions and issues, ranging from processes within the reservoir itself (flow, mineral precipitation, channelling, etc.) to the access and energy systems (wells, pumps, heat exchangers, surface tubing insulation, etc.).

	- a. Well integrity
		- i. Induced primary loop corrosion (e.g., acidic waters corroding steel), both internal and external to the well casing, potentially exacerbated by high temperatures and electrochemical corrosion
		- ii. Primary loop internal mineral precipitation—scaling (increasing pipe friction losses)
	- b. Reservoir condition
		- i. Mineral and rock dissolution (e.g., dissolving of gypsum) or other alteration (e.g., induced clay mineral swelling from changing geochemistry).
		- ii. Microbiologically induced pore blockage or corrosion through the generation of biofilms or weak organic acids that act on minerals (mainly carbonates)
		- iii. Chemical changes and processes induced by flow, ΔT, and ΔP, leading to solubility gradients that trigger dissolution and precipitation
		- iv. Flushing of pipe corrosion and precipitated mineral particles into the reservoir, generating blockage of flow paths
	- a. Rock mechanics issues: stress-strain (σ–ε) processes and their effects on porosity and permeability, particularly in systems dominated by fracture flow and susceptible to thermoelastic impacts
	- b. The transition from laminar to turbulent flow regime near the well at high flow rates (short-term impacts on production and injection phases)

The long-term result of disequilibrium processes in the reservoir may be evidenced as pore throat blockage and fracture aperture reduction, leading to reservoir fluid conductivity changes, including permeability reduction, channeling of flow, and related processes. If these changes are properly understood and predicted, project planning becomes easier, allowing answers to questions related to the drilling of new wells, workover scheduling for impaired wells, adding more injection or production wells, developing a new productive horizon (deeper or shallower), changing the electricity/heat output ratio, and so on.

To optimize energy production, increase the project life-time, enhance the fiscal outcomes and minimize environmental impacts, many issues associated with heat extraction from a geothermal well-doublet system need to be considered; these issues are addressed in this study.

#### **2. Chemical Condition**

Questions arise as to whether deep geothermal systems can be preserved while retaining their production quality, achieving the desired E over the project life. We note ˙ that the best outcome is keeping E constant, but a project can also be designed on the basis ˙ of a gradually declining E, as long as it can be realistically predicted. Answers depend ˙ on water-rock reactions and reservoir condition changes over the project life, as chemical reactions may cause flow-path changes (precipitation in flow paths) and potentially alter the heat transfer efficiency between rock and water. To predict changes in permeability or fracture conductivity, chemical effects need to be investigated for the project's expected life-span [29].

#### *2.1. Well Integrity*

Hot (>60 ◦C) geothermal fluids are often associated with dissolved carbon dioxide (CO2) and hydrogen sulfide (H2S), which makes the saline brine slightly acidic, providing an aggressive thermochemistry environment [30], leading to severe internal corrosion of steel goods (casing, pipes, heat exchangers, etc.) through the sulfide reaction pathway, shown in Figure 2. Note that well integrity discussions include physical damage (shear, cracking) as well as applied mitigation measures to control corrosion. Corrosion/scaling damage reduces the system efficiency by changing the heat production capacity (Q) and adding corrosion products to the circulation system, potentially reducing reservoir permeability. Well integrity will be affected, shortening productive life [29] and generating various effects leading to efficiency loss, increased maintenance costs, and surface damage [31]. At a cost, anti-corrosion materials can be added to the system, and the choice of material depends on ΔT, ΔP, and formation water chemistry.

**Figure 2.** Schematic illustration of the process of iron dissolution and sulfide precipitation in a CO2/H2S aqueous system (modified from [29]).

> External steel casing corrosion and cement degradation may occur in an accelerated manner if the formation temperature is increased (Figure 3). Existing outside-the-casing groundwater chemistry in the upper geological section will change from chloride-based to sulfate- or bicarbonate-based as the surface is approached in most sedimentary basins. Heating the cemented steel casing system sets up a large electrolytic cell where corrosion is accelerated by the temperature increase from rising hot fluids. External pitting leads to casing perforation and integrity loss.

**Figure 3.** External steel casing corrosion from acidic water attack [32].

In certain chemical conditions, specifically, the combined presence of archaebacteria, CaSO4·2H2O, and a slow source of CH4 in a stratum, acidic conditions (H2S generated) are created and external steel corrosion accelerated by both the +ΔT and the increasing acidity. We note that the methane may even be traveling up the geothermal production well outside the casing [33].

#### *2.2. Reservoir Condition*

Mineral dissolution and precipitation, ΔT and ΔP, and steel goods corrosion are the main issues related to chemical disequilibrium in deep geothermal systems [34]. Geothermal fluids have diverse chemical components due to the diversity of geological settings (e.g., silica-rich, carbonate-rich, chloride brines, etc.). The recharging water source and associated dissolved gases also define the chemical characteristics of the geothermal fluid, and these may change with time because of the perturbation caused by production and injection. The precipitation of minerals or induced corrosion/scaling leads to small particles carried in suspension into the porous reservoir, altering the system's porosity and permeability. Bächler and Kohl [35] investigated such changes in a deep geothermal system using coupled Thermal-Hydraulic-Chemical (THC) modeling. The study found that the system's efficiency decreased more rapidly with more mineral dissolution at the start of the project, and it slowed down over a longer time. These results were attributed to the disequilibrium of the system at the start, then changes in the equilibrium condition over time [35]. For production and service life calculations, flow pathway impairment from particles and precipitates is an important disequilibrium effect.

Microbially induced (mediated) corrosion, MIC [36], and fouling affect geothermal energy projects (as well as other industrial activities); designing and implementing protection requires multiple knowledge perspectives—material science and metallurgy, electrochemistry, biochemistry, and microbiology [37]. Salas et al. [38] showed MIC's impact on the deterioration of geothermal power plant structures (e.g., vapor ducts and cooling tower supports), noting bacterial activity at temperatures as high as 140 ◦C. A MIC problem in the primary surface loop may migrate to the subsurface as corrosion products, and temperature changes impact the reservoir's chemical conditions and change the system's flow properties [39]. Tubular goods may experience general (e.g., steel thinning) or localized corrosion (e.g., pitting) [40]. Uniform corrosion [41], pitting, crevice, and contact (galvanic) corrosion [42–44], microbially mediated corrosion [45], and oxygen corrosion [38] are all common in geothermal systems.

Thermohydrochemical reservoir processes are impacted by stimulation efforts (e.g., hydraulic fracturing with gels, acid injection, etc.), and injection/production leads to ΔT and ΔP, triggering mineral precipitation. Even reservoir rock and fluid heat conductivities can be altered by ΔT and ΔP. Blöcher et al. [46] showed that significant changes occur in the porosity and permeability of a sandstone reservoir when the effective stress (σ) is increased, as shown in Figure 4.

**Figure 4.** Decrease in measured porosity and permeability of Flechtinger Sandstone geothermal reservoir because of increased effective stress (here called effective pressure) [46].

Remedial/preventive methods include cleaning, jetting, and waste removal in well workovers, and soft acidizing during stimulation to remove carbonate precipitates, with remediation processes guided by leak-off and injection step-rate testing and data from other chemical and physical monitoring methods [47–49]. Note that soft acidizing of geothermal systems is taken to mean injecting the same acid volume as used in conventional oil and gas acid treatments but at a much lower rate [48,50]. Salimzadeh and Nick [51] developed a two-way chemical coupling method to include mineral dissolution in Thermal-Hydraulic-Mechanical (THM) coupling, generating a coupled Thermal-Hydraulic-Mechanical-Chemical (THMC) approach to explain noticeable changes in reservoir properties arising from chemical disequilibrium.

#### **3. Geomechanical Conditions**

Fluid injection into a geothermal well system can induce stress changes at a scale through hydromechanical, poroelastic, and thermoelastic mechanisms. In this section, we evaluated the geomechanics responses for cold water injection into a geothermal double-well system. Following geophysics symbology, we use "S" to mean total stress, which at a point is the sum of the effective stress σ and the pore fluid pressure, viz.: Sv = σ<sup>v</sup> + Pp; SHmax = σHmax + Pp; Shmin = σhmin + Pp. Injecting cold water into hot rock perturbs the stress in two ways: through ΔP and temperature ΔT. As a preliminary example, assuming a uniform planar and linearly elastic stratum, the effective horizontal stress change (Δσ<sup>h</sup> − MPa) for a uniform reservoir wide temperature change (ΔTr − K) is

$$
\Delta\sigma\_{\rm h} \approx \frac{\Delta\Upsilon \times \to \times \alpha\_{\Gamma}}{1 - \upsilon} \tag{2}
$$

where the elastic properties defining the volume change (ΔV) are Young's modulus (E − GPa), Poisson's ratio (ν), and the coefficient of thermal expansion (α<sup>T</sup> − <sup>K</sup>−1); for a stiff rock E ≈ 60 GPa, <sup>α</sup><sup>T</sup> ≈ <sup>10</sup>−<sup>5</sup> <sup>K</sup>−1, and <sup>ν</sup> ≈ 0.2, so that a <sup>Δ</sup>T of −30 K can lead to a large effective stress change, Δσ<sup>h</sup> ≈ −22 MPa in this example. In cases with steep temperature gradients from cooling (across boundaries of large permeability contrast, for example), stress changes can readily induce shear slip of pre-existing fractures or small faults and exceptionally trigger the shear yield or tensile rupture of intact rock. This is also impacted by pressure changes, and to evaluate shear yield, it is necessary to introduce a shear slip criterion.

The potential for fracture slip or shear yield of intact rock from combined ΔT and ΔP effects from the fluid injection can be evaluated using simplified (i.e., linearized) Mohr-Coulomb (MC) yield criteria [52]

$$
\pi = \mathbb{C} + \sigma\_n \tan \varphi \tag{3}
$$

where, τ is the shear stress at yield (shear slip) along a plane, usually, a rock weakness plane in naturally fractured rock masses, and σ<sup>n</sup> is the normal effective stress acting across the slip plane. C and ϕ are the cohesion and internal friction angle, respectively; that is, the shear strength parameters of the slip surface (e.g., fault or fracture surface). To apply the MC criterion, the in-situ effective stress tensor, the slip plane shear strength, and the slip plane orientation are needed. Determining the effective stress tensor (σij) at the region of potential slip is a coupled THM problem.

Fluid injection changes stress through hydromechanical (ΔP − Δσ), poroelastic (ΔP − ΔS), and thermoelastic (ΔT − Δσ) volume changes that generate three-dimensional strains (εij) and therefore stress changes (ΔSij and Δσij). Of course, Δσij, in turn, leads to more ΔV; hence, the term "coupled" is used: no process is independent of the others. The first mechanism accounts for the bulk compressibility effect ΔV, the second for the porous rock ΔV, and the third for the thermally induced ΔV. Quantification of ΔV requires coupling heat conduction and convection with fluid flow and stress-strain analysis. In general, in the low permeability host rock, higher pressures reduce the frictional strength (Equation (3)) of a slip plane by decreasing the effective normal stress (σij = Sij − Pδij), facilitating yield. In a Mohr stress diagram (Figure 5a) using total stresses, a pressure increase shifts the

yield envelope (black line) into the Mohr circle (the stress state), and yield is said to have occurred when stresses lie outside the criterion limit.

**Figure 5.** (**a**) 3D Mohr circles showing stress condition and slip tendency due to fluid injection in the deep geothermal reservoir Groß Schönebeck in the North German Basin. The red line is the Coulomb failure criterion when pore pressure has been elevated by 7 MPa, and (**b**) The lower hemisphere stereonet plot illustrates the slip-tendency (ratio of resolved shear to normal stress) and stress direction.

Because stiff, brittle earth materials store elastic strain energy but are strain-weakening, a sudden stress drop (weakening) accompanies shear slip, and a very small fraction of the released strain energy is evidenced as strain waves (seismic waves). Although the fluid-injection fracture-slip-induced deformation is not high, generally far less than a centimeter, Δσij may also trigger and reactivate nearby faults that are near critical stress conditions, causing felt seismicity and casing shear problems [53]. Injection-induced earthquakes are usually small with moment magnitude < 0. However, there are cases where large magnitude (>1) injection-induced seismicity has been recorded [54]. Two groups of field parameters can impact the magnitude and rate of injection-induced seismicity: controllable operational parameters that include fluid injection pressure, rate, temperature, and volume [54]. Uncontrollable subsurface parameters include the initial state of stress and pore pressure, size and density of pre-existing faults/fractures, fault/fracture orientation and shear strength, and other geomechanics parameters [55].

Figure 5 shows the effect of the pressure increase on the fault reactivation in the deep doublet geothermal Groß Schönebeck project in the North German Basin, where fluid-injection-induced seismicity is recorded [56]. The stress state at the injection depth of 4035 m is a normal faulting regime with Sv = 100 MPa, SHmax = 98 MPa, Shmin = 55 MPa, and Pp = 43 MPa [56]. Faults striking in the SHmax direction with dip angles of 45–75◦ are critically stressed; increasing pore pressure by less than 7 MPa increases the risk of induced seismicity on two fault sets: the NW-SE set (rectangle), and the E-W set (Mohr circle). Shear slip of natural fractures in strong rocks is usually a dilatant process, so the rock mass permeability is enhanced in certain directions, which in turn will alter the fluid flow patterns, and probably increase the permeability anisotropy in a rock mass.

The second mechanism, poroelasticity, governs fluid pressure change effects on insitu states of stress and rock mass deformation [52,57]. Geothermal SedHeat wells are drilled into porous, permeable reservoirs, and ΔP causes volume changes that lead to stress changes, altering porosity and permeability. Engelder and Fisher [58] demonstrated that stress/pore pressure changes in a normal faulting regime as follows

$$\frac{\Delta \mathbf{S}\_{\text{hmin}}}{\Delta \mathbf{P}} = \alpha \frac{1 - 2\upsilon}{1 - \upsilon} \tag{4}$$

where α is Biot's poroelastic coefficient, for this relationship, it is assumed that vertical stress remains constant.

In this study, we evaluated the pore pressure stress coupling at different tectonic stress settings employing 3D geomechanics modeling. According to Anderson's classification scheme [52], the three states of stress are normal (Sv > SHmax > Shmin), strike-slip (SHmax > Sv > Shmin), and reverse faulting regimes (SHmax > Shmin > Sv), where S values are the principal compressive stresses. One-way hydromechanical coupling is performed to evaluate the fluid injection effects on the stress tensor using the finite element platform Visage™ Geomechanics Simulator. We considered a hypothetical reservoir located at a depth of 2250 m where Pp is hydrostatic (≈11 MPa/km), and Sv increases at a rate of ≈26 MPa/km. Simulation boundaries were sufficiently distant to have no effects on the deformations in the affected (reservoir) zone. This model was subjected to a +ΔP of 8 MPa.

Figure 6 shows the results for (a) normal faulting, (b) strike-slip, and (c) reverse faulting stress regimes, respectively. Pre- and post-injection stress tensors at the point of fluid injection are illustrated by 3D Mohr circles. In all cases, injecting fluid shifts the Mohr circles (stress magnitudes) leftward toward the yield criterion, indicating that the stress state is closer to the yield criterion. Note that in the case of the normal fault, incipient slip is apparent for the values chosen (Pp = 35 MPa). In the normal faulting regime, the effective maximum and minimum principal stress difference (σ<sup>v</sup> − σhmin) decreases as Pp is increased, but not for the other two regimes; notably, in different states of stress, the local stress tensor components respond differently to fluid injection.

**Figure 6.** Mohr's Circle in three dimensions showing the effect of fluid injection (for pressure difference, ΔP = 8 MPa) at different states of stress: (**a**) normal faulting, (**b**) strike-slip, and (**c**) reverse faulting stress regimes.

The third ΔV mechanism, thermoelasticity, leads to slow rock shrinkage with cold water injection [52], inducing stress changes. The thermally induced stress magnitude (ΔσT) effect was illustrated above (Equation (2)). We re-iterate that the stress redistribution process is driven by volume changes governed by ΔT and the host rock elastic properties [59]: ΔV=V·ΔT·αT, where the bulk rock thermal-expansion coefficient is αT. As shown above, in stiff rocks (high E), Δσ<sup>T</sup> effects can be very large. Cooling occurs mostly convectively in permeable rocks, but in bounding low-permeability rocks, only conduction obtains. Therefore, strong local time-dependent differences in thermal expansion can take place, leading to shear stress concentrations at bounding interfaces (usually bedding), leading to shearing and induced seismicity. Figure 7 demonstrates how cool fluid injection changes the local state of stress around the cooled zone. Thermally-induced contraction decreases σij in the cooled zone, and redistributes stresses around to maintain mechanical stress equilibrium. A similar (poroelastic) effect is expected to take place with ΔP, but the two processes, ΔT and ΔP, have widely different characteristic time scales. This implies that the rock deformation and stress changes are likely to be dominated by ΔP in the early time of a geothermal project, whereas ΔT effects become increasingly more important during the late time [60].

**Figure 7.** Induced thermal stress resulting from cold fluid injection into a geothermal reservoir. The color scale is the relative temperature range from 60 to 100 ◦C.

Hence, the zone being cooled relative to the rock mass loses compressive stress from thermoelastic shrinkage; this lost stress is transferred as increased stress to the surrounding, uncooled rock to maintain stress equilibrium (Figure 7). This stress redistribution generally increases local shear stresses, but the zone of maximum shear stress increase is distant from the injection well, near the region of the thermal front, which can be quite sharp because

the reservoir rock shrinks quickly (convectively), whereas the overburden shrinks very slowly (conductively).

In geothermal well systems, increasing the flow rate of the geothermal fluid can lead to the generation of a local turbulent flow regime, particularly around the injection well due to pumping the fluid flow with high pressure, affecting efficiency. In the short-term, the high flow rate of the injection fluid can lead to an increase in thermal energy production, although it has the inverse effect in the long-term. In addition, the turbulent flow near the injection well can cause more stress and, consequently, damage the well because of the induced stress created over time.

#### **4. Heat-and Flow-Transport Challenges**

#### *4.1. System Thermal Issues*

Eventually, loss of capacity to produce sufficient fluid rates at high enough temperatures to sustain commercial operation defines a geothermal reservoir's life-time [46]. The reservoir characteristics of porosity, permeability, initial temperature, and bulk heat capacity are performance indicators of a geothermal injection/production well system. These indicators dictate daily energy production and reservoir life-time. The recovery factor is the amount of producible energy compared to the total available energy, and it depends on the use of the heat at the surface (electrical power, power + heat, seasonal heat storage, etc.). Pure power production leads to a lower energy recovery factor. However, if the geothermal system is cooling, useful heat can likely be produced long after the period of electricity generation. Capital cost, operational cost, and the geothermal system performance over time define commercial viability [47]. The performance of a geothermal well system (commercially viable life-time and energy recovery) is related to both the natural physical parameters listed above and human-controlled factors. Transitioning a decaying geothermal system from power generation to sensible heat use makes great economic sense because the capital investment payback is dominated by the early years, when electricity is being generated. For long-term sustainability, this means there must be additional SedHeat systems added on episodically as the initial systems decay because of the slow drop of heat energy provision.

Human controlled parameters are [48]: (1) discharge rate from the reservoir, (2) injection rate and fluid input temperature, and (3) well spacing. Production, injection, and well spacing effects have been studied by Kazemi et al. [8] to show their importance on the lifetime of a doublet geothermal system. Porosity, permeability, chemical components of pore fluid, initial temperature, reservoir minerals, and their chemical components and mechanical parameters, reservoir thickness, and the presence of shale layers and their thicknesses are the main natural subsurface physical parameters [48,49]. Thermal depletion happens by continuing injection of cooler water until the production well shows thermal breakthrough and decline in temperature, impacting the efficiency of the doublet geothermal system. Various means of extending life-span and project utility might include episodic electrical power production through integration of other energy sources (wind, solar, diesel, etc.), developing such a large well spacing that life-span is increased (but recovery efficiency may drop), and by thermal recharging from industrial waste heat sources or by use of excess solar and wind energy when available. Solar thermal collectors can be installed instead of photovoltaic systems, and these can feed the hot fluids into the system before the heat exchanger or afterward [61,62], depending on the time of year, the electrical power demand, and other factors. In all cases, the thermal capacity retention time of a project depends significantly on human-controlled parameters, given a set of physical system parameters [50].

Production-recovery cycles comprise an important process in any geothermal well system and need to be considered in any life-time assessment and economical assessment of such a system. A case study on a Riehen reservoir exploited by a doublet geothermal system showed that the recovery of the reservoir is closely related to the cyclic production-recovery plan [50], suggesting that shorter cycles allowed more thermal energy to be produced. If there is some component of natural fluids recharge (a leaky reservoir), beneficially encouraging such recharge involved not only rate management but pressure management. Thermal breakthrough will occur sooner if production continues without letting the reservoir regain some heat through conduction from the surrounding rocks. The exploitation of the doublet aquifer is also affected by the pressure gradient in the aquifer, changing heat retention in the reservoir [63]. Table 1 summarizes some parameters for more effective heat production from a doublet geothermal system.

**Table 1.** Effective parameters for heat production from a doublet geothermal system.


The main thermodynamic/commercial objective of a geothermal system must be economically viable heat recovery efficiency, defined as the ultimate recovered heat over the reserved heat in the aquifer [64–68], subject to commercial criteria. In order to achieve optimum efficiency from a geothermal system, the main goal should be to optimize heat recovery from the reservoir to service both the electrical grid needs and to provide useful heat for habitat and industrial needs. Most geothermal power plants are designed to provide electrical energy/heat to the surrounding area, which often is a local detached grid [69–72]. If geothermal electrical power can be used in a large grid to provide special services such as peak shaving or load displacement, the value of the energy provision increases.

#### *4.2. Reservoir Geological Aspects*

A first-order factor in the productivity of a geothermal SedHeat well array is the stratigraphy and flow properties of the reservoir. Figure 8 is a representation of a doublet in a typical SedHeat situation.

The target stratum is bounded by impermeable strata: U shale is the upper advective flow boundary, L shale is the lower advective flow boundary. Given typical permeability contrasts in geological media, it is usually assumed that the SedHeat fluid reservoir is bounded by impermeable strata, such that advective flow is confined to the more permeable intervals. In litharenite, SedHeat reservoirs (sand-silt-shale), the productive interval contains granular beds (sandstones) of differing permeability; the most common distribution is to find higher values in the lower beds, grading upward to lower permeability in the uppermost beds (also characteristic of transgressive sedimentary regimes).

Commonly, there are intraformational low-permeability strata, deposited under more quiescent flow conditions, called silts or shales; they are of much lower permeability than the low k sandstone, usually by more than two orders of magnitude. Furthermore, because of processes in the original depositional environment, these low-permeability strata may be continuous at the doublet spacing scale (shale C), or they may be laterally discontinuous (shale A and shale B). The complete disposition of the various strata between the doublet pair is difficult to specify, even with excellent seismic investigation; usually, stratigraphic details are known only from the geophysical logs available from the production and injection wells.

**Figure 8.** Geological disposition factors in doublet efficacy.

One may assume that in the stratigraphic disposition shown, the permeability of the lowermost sandstone is several times that of the central sandstone, which in turn is several times that of the upper sandstone (a permeability range for the sandstones might be from 1 D to 0.1 D). The silty shale interbeds (labeled "shale") have permeability values less than 0.001 D, so they may be considered functionally impermeable in the context of advective fluid flux. Thus, heat in these strata (water and mineral phases) can migrate to the sandstone only by conductive heat transport, and is thence taken to the production well advectively (blue arrows). Similarly, the heat in the bounding strata (U-, L-shale) can only flow conductively and be exploited when the advectively flowing fluid in the sandstone drops below the original temperature.

Another issue in a natural, unbounded reservoir is that flow is not constrained by lateral no-flow boundaries. Assuming an original pressure Po and similar pressure build-up and draw-down cones around each well (blue line), part of the injected fluid will dissipate into the farfield to the right, driven by the positive ΔP, and at the production well, part of the warm fluids will come from the farfield to the left, driven by the pressure sink -ΔP. The heterogeneity that is characteristic of all sedimentary rock geothermal reservoirs invariably decreases the energy recovery and the energy recovery rate. The most optimistic recovery scenario is the isotropic, constant permeability case, with the further assumption that the local permeability is not impaired by geochemical processes or fine-grained solids migration into the reservoir.

Finally, we point out that in geothermal reservoirs that are dominated by flow in natural fractures, there is an important short-circuiting mechanism associated with thermal effects: maximum fracture aperture increase occurs in the fracture with the most flow, reinforcing the flow rate [73]. It remains unclear if this process can be somewhat suppressed by flow management and by thermal recharge [74]. Assessing the commercial viability of a geothermal SedHeat reservoir that involves the injection of cold water is and will remain a challenging task. Currently, SedHeat projects in laterally extensive reservoirs that do not require re-injection of cold water seem commercially viable because the fluid entry temperature remains constant from continuous recharge. Injection/production projects based on well-doublets or variants require more careful evaluation.

#### **5. Conclusions**

Geothermal energy appears to be a green, generally available, potentially constant, and stable source of energy. If economically viable, it is a practical and sustainable solution for supplying energy needs (power or thermal) and can decrease environmental challenges and threats (e.g., air pollution, global warming, climate change). Geothermal well systems can be used to generate both electrical and thermal energy. Despite the advantages listed here, these systems will evidence challenges over time, including heat conductivity, chemical, and mechanical issues, which will affect their operational efficiency, as well as costs. Based on the literature and studies conducted, the following results can be concluded:


Knowing the issues discussed in this article will help designers/engineers to better design and implement geothermal SedHeat systems, based on the local geological characteristics of the resources and geographical and climate conditions, resulting in improving operations and enhancing efficiency of the systems as well as reducing their associated costs and risks.

**Author Contributions:** All authors contributed to the writing and editing of the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study received no external funding.

**Acknowledgments:** The authors would like to thank Schlumberger for making the Petrel E&P Platform available for this study.

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

#### **Nomenclature**


#### **Abbreviations**


#### **References**

