Numerical Investigation on the Inﬂuence of Surface Flow Direction on the Heat Transfer Characteristics in a Granite Single Fracture

: As ﬂuid passes through the fracture of an enhanced geothermal system, the ﬂow direction exhibits distinct angular relationships with the geometric proﬁle 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 inﬂuence of the ﬂuid ﬂow direction on the heat transfer characteristics of granite single fractures and the accuracy of the numerical model was veriﬁed by experiments. Results demonstrate a strong correlation between the distribution of the local heat transfer coefﬁcient and the fracture morphology. A change in the ﬂow direction is likely to alter the transfer coefﬁcient value and does not affect the distribution characteristics along the ﬂow path. Increasing injection ﬂow rate has an enhanced effect. Although the heat transfer capacity in the fractured increases with the ﬂow rate, a sharp decline in the heat extraction rate and the total heat transfer coefﬁcient is also observed. Furthermore, the model with the smooth fracture surface in the ﬂow direction exhibits a higher heat transfer capacity compared to that of the fracture model with varying roughness. This is attributed to the presence of ﬂuid deﬂection and dominant channels.


Introduction
Geothermal energy is a clean and environmentally friendly renewable energy source with a wide distribution range, large reserves, and long duration [1][2][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 CO 2 ) [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].
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.

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). 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.

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). 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 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.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 4 of 20 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.  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: Momentum equation: Energy equation: where (kg/m 3 ) and (kg/m 3 ) 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; (W/(m·°C)) is the fluid thermal conductivity; (Pa s) is the fluid viscosity.  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: Momentum equation: Energy equation: where ρ w (kg/m 3 ) and ρ r (kg/m 3 ) 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; λ w (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: Specific heat capacity of water: Thermal conductivity of water: Water density: 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.

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. 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: 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.

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.  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  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.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 6 of 20 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.

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. Calculations were performed following the unsteady double precision method and the second-order implicit mode to improve calculation accuracy. The SIMPLEC (semi-implicit 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 −5 and the time step size for the transient simulations time was equal to 1 s. The total solution time was taken as 30 min.

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

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.

Model Case (α =)
Rock Temperature, T 0 ( • C) Flow Rate, V 0 (mL/min) 10 20 30 40 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 −5 and the time step size for the transient simulations time was equal to 1 s. The total solution time was taken as 30 min.

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. 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. 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.  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. 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. 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.
where Aout (m 2 ) 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.  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 T outlet at the fracture outlet for a given time was calculated as follows: where A out (m 2 ) 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.
where Aout (m 2 ) 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.

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]: where M (kg/s) is the flow rate; H inlet (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.

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]: where q f (W/m 2 ) is the average heat flux on the inner wall of the fracture; T f ( • C) is the average temperature on the inner wall of the fracture; T w ( • C) is the average temperature of the water in the fracture. These parameters can be calculated from the numerical simulation results as follows: where A f (m 2 ) is the surface area of the fracture; V f (m 3 ) is the volume of the fracture.

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: where q f x (W/m 2 ), T f x ( • C), and T wx ( • 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: 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: where d e (m) is the equivalent hydraulic diameter of a channel with arbitrary shape. 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 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  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 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.

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

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 0 = 90 • C compared to that with initial temperature T 0 = 70 • C. The maximum Nusselt number at α = 90 • is 15.8 (Figure 13b).

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 T0 = 90 °C compared to that with initial temperature T0 = 70 °C. The maximum Nusselt number at α = 90° is 15.8 (Figure 13b). 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: 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. 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: Figure 14 presents the cumulative heat extraction rate Q sum (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 (∆Q sum , 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. 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 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 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 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 Figure 15. Variation of heat extraction rate with time at several initial flow rates and initial temperatures. 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 • . 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 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  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. 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.  z-axis. Similar to the fracture profile, the distribution curve of the local heat transfer co ficient is characterized by peaks and troughs. The local heat transfer coefficient curve α = 0°, α = 30°, and α = 60° reveal that the local heat transfer coefficient at the peaks of fracture profile is consistently greater than that at the troughs. Moreover, increasing injection flow rate not only improves the local heat transfer coefficient, but also has enhancing effect that makes the difference between the peaks and troughs of the local h transfer coefficient more evident. The enhanced effect can be observed via the increa slope of the local heat transfer coefficient curve. Figure 19 shows the enhancement eff of increasing the flow rate on the convective heat transfer intensity.    Figure 20 presents the distribution of the local heat transfer coefficient, wall temp ature, and average flow rate on the fracture profile at z = 0.05 m. The curves of the avera flow rate and local heat exchange coefficient become increasingly complex as α increas from 0° to 90°. Since the fracture model is constructed by translating along the y-axis, t distribution of fracture aperture in each direction of the model is not uniform. Therefo a position with a higher absolute gradient results in a smaller fracture aperture, whi results in a low flow rate as the water flows through the path of the small aperture. contrast, fractures with larger aperture will form a dominant channel where the flow ra is greater. Compared with the temperature, the difference of flow rates at fractures w different apertures is the main factor causing uneven distribution of local heat trans coefficient. In addition, a strong positive correlation is observed between the local he exchange coefficient and the flow rate distribution, which indicates that the flow rate the main factor affecting the heat transfer intensity.   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 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.

Local Heat Transfer Coefficient
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.

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.

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 in-

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.

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 Q sum 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.