A Novel Procedure for Coupled Simulation of Thermal and Fluid Flow Models for Rough ‐ Walled Rock Fractures

: An enhanced geothermal system (EGS) proposed on the basis of hot dry rock mining tech ‐ nology has become a focus of geothermal research. A novel procedure for coupled simulation of thermal and fluid flow models (NPCTF) is derived to model heat flow and thermal energy absorp ‐ tion characteristics in rough ‐ walled rock fractures. The perturbation method is used to calculate the pressure and flow rate in connected wedge ‐ shaped cells at pore ‐ scale, and an approximate analyti ‐ cal solution of temperature distribution in wedge ‐ shaped cells is obtained, which assumes an iden ‐ tical temperature between the fluid and fracture wall. The proposed method is verified in Barton and Choubey (1985) fracture profiles. The maximum deviation of temperature distribution between the proposed method and heat flow simulation is 13.2% and flow transmissivity is 1.2%, which in ‐ dicates the results from the proposed method are in close agreement with those obtained from sim ‐ ulations. By applying the proposed NPCTF to real rock fractures obtained by a 3D stereotopometric scanning system, its performance was tested against heat flow simulations from a COMSOL code. The mean discrepancy between them is 1.51% for all cases of fracture profiles, meaning that the new model can be applicable for fractures with different fracture roughness. Performance analysis shows small fracture aperture increases the deviation of NPCTF, but this decreases for a large aperture fracture. The accuracy of the NPCTF is not sensitive to the size of the mesh.


Introduction
Geothermal energy, which is clean, renewable, and widespread, has been developed and utilized in many countries, and considerable attention has been given to exploiting thermal energy from hot dry rock (HDR) 3-10 km underground, using the enhanced geothermal system (EGS) to exchange and loop heat in/out of effectively connected systems of underground fractures [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. For example, the Landau plant in Germany [15], the Soultz plant in France [16], and the recently commissioned Geodynamics Habanero pilot plant [17] demonstrate the feasibility of EGS running. Figure 1 shows the traditional heat extraction process from the geothermal system, where the vast rock fracture system is the main path for heat fluid flow. However, prediction of EGS service life and thermal energy extraction encounters considerable challenges due to limited knowledge of the heat exchange behavior between flowing fluid and high temperature rock [18][19][20][21]. Natural rock fractures are characterized by irregular shape, rough surface, and even contact asperities, which play important roles in the hydraulic and heat transfer properties of fractures. Because of the rough fracture surface, there may be a 1-2 order of magnitude error of permeability calculated by the widely used cubic law, which assumes the fracture is a smooth parallel plate [22][23][24]. Zou et al. [25] and Wang et al. [26] used a 3D Lattice Boltzmann method to predict the fluid flow behavior in fractures and found that the secondary roughness enhanced the local complexity of velocity distribution and led to nonlinear flow behavior. Similar studies were also reported by Yin et al. [27] and Xiong et al. [24]. In addition, the rough fracture surface also causes a decrease in heat transport intensity. Recently, many efforts [28][29][30] paid attention to the study of the effects of fracture roughness on heat transport behaviors in fractures, based on laboratory experiments and numerical simulations. For example, He et al. [28] found that heat transport characteristics mainly depended on fracture surface roughness, followed by aperture and flow rate, combining the experimental and numerical modeling approaches. Huang et al. [29] conducted seepage and convective heat transport experiments, and demonstrated that large roughness in the direction perpendicular to flow would increase the capacity of heat transport. Ma et al. [30] adopted numerical and experimental approaches to improve understanding of the heat transfer characteristics of water flowing through rough fractures, and results showed the total heat transfer coefficient increased with the increase in the value of the joint roughness coefficient (JRC). However, there are no analytical models describing the heat transport properties between fluid and high temperature rock, considering fracture surface geometry. In previous research studies, the parallel plate is traditionally used to simplify the real rough fracture surface, and many models can be observed in the literature, such as the Gehlin and Hellstrom model [2], Gringarten et al. model [31], Cheng et al. model [32], Martinez et al. model [33], Zhao model [34], Yost and Einstein model [35], and Yan and Jiao model [7] which have inevitable relative error comparing with the real heat transport process in rough-walled fractures.
In view of analytical heat transport models for rough fractures being rare, this paper aims at developing a proposed novel procedure for coupled simulation of thermal and fluid flow (NPCTF) models to understand the heat transport process in rough fractures, based on representing the measured fracture geometry with a series of connecting wedges formed using adjacent apertures in longitude direction along the fracture plane, which is a common approximation approach for fracture geometry [36]. The performance of the proposed NPCTF is assessed by comparing predictions with results from numerically solving the Navier-Stokes and energy equations, using COMSOL code in rock fracture profiles obtained by a 3D stereotopometric scanning system, which have different relative surface roughness and aperture.

Fluid Flow Model
Flow behavior in natural rock fractures has been studied through considering fracture walls as parallel plates, saw-tooth shaped walls, and sinusoidally varying walls [30]. Among these approximates, fractures are found to be well described by a series of connected wedges. Wang et al. [37,38] used a wedge-shaped fracture to obtain the perturbation solution of pressure ∆p and flowrate Q under the pressure boundary condition. The perturbation parameter ε is selected as the relative aperture variation along the wedge length l in his study. This section provides a brief description of the approach rather than a detailed derivation. A more detailed perturbation derivation can be viewed in Wang et al. [37,38]. By perturbation analysis, the pressure difference can be made dimensionless and expressed as expanded series with a small parameter ε: where X is defined as X = εx/bm; ω is the dimensionless absolute aperture variation defined as ω = |a|/bm, where a is the difference between the upper and lower wedge edge; B is the dimensionless aperture defined as B = c/cm, c is the fracture local aperture, and cm is the mean aperture; B′ is the first derivatives of B with respect to X; R is the Reynolds number defined as R = ρQ/μ; Q is flow rate per unit width of fracture and μ is fluid viscosity; ∆P is the dimensionless pressure difference given as ∆P = ∆p/∆pm and ∆pm is the pressure difference of flow through a fracture with a uniform aperture defined from the cubic law as: Noted that only the first two order terms in Equation (1) are used which is enough accuracy for nonlinear flow simulation. Because the quadratic polynomial form, called the Forchheimer law, can well describe nonlinear flow behavior in a single fracture [39,40], the highest order term is the second order term.

Heat Transport Model
The steady heat fluid transport differential equation in the fracture, including advection, conduction, and convection terms for the fracture walls, can be expressed as [34]: where ρw is the fluid density; cw is the specific heat of fluid; Kw is the fluid thermal conductivity; Kr is the rock thermal conductivity; v is the steady flow velocity. b is half aperture of the fracture; Tf is the bulk temperature of the fluid; Tr is the temperature of the reservoir rock matrix. In this equation, the temperature at the rock fracture surface is identified as that of the bulk fluid. The steady heat conduction in the rock is governed by [32,34]: which assumes that the heat conduction is two dimensional, perpendicular to the fracture plane.
The boundary conditions associated with the heat flow in the fracture are: where R is the height of the rock. The identical temperature between the fluid and fracture surfaces is expressed as: , , where ft(x) and fl(x) are the upper fracture wall and lower fracture wall, respectively. Zhao et al. [34] used this assumption to obtain a solution of temperature in the plane fracture. This solution has good agreement with the experiment results. However, there are not any temperature solutions for the fracture wedge, which is the basic element in rough fractures. Therefore, a solution for the fracture wedge would be developed based on the above heat fluid transport equations.
(1) Symmetric fracture wedge For the symmetric fracture wedge (Figure 2a), the wedge wall (upper fracture ft(x), lower fracture fl(x)) can be described as: The aperture b can be determined as ft(x) − fl(x) according to Equations (12) and (13). If thermal diffusion in the fluid is not considered, the boundary conditions in Equations (10) and (11) are substituted into Equations (5) and (6); a temperature distribution equation of fluid can be expressed as: The solution to Equation (14) considering the temperature boundary in Equations (7)-(9) is: (2) Asymmetric fracture wedge For the asymmetric fracture wedge (Figure 2b), the wedge wall can be described as: Similarly, the solutions to Equations (5) and (6) for the asymmetric wedge with boundaries in Equations (16) and (17) The temperature mathematic solutions for the fracture wedge were first proposed in the literature based on a few assumptions which have well-defined parameters.

Model Description
As presented above, the fluid flow model (Equation (1)) and heat transport model (Equations (18)-(20)) were obtained in the fracture wedge. Each fracture wedge was connected to another and constituted a rough-walled fracture ( Figure 3a). According to the flowrate equilibrium principle and fluid flow model (Equation (1)), the flowrate and pressure distribution in the entire fracture can be determined (seen in Figure 3b). If the flowrate distribution is known, the temperature solution in each connected fracture wedge would be obtained using the heat transport model (as Figure 3c). When the pressure and temperature condition changes, the density of water (ρw) is no longer a constant value; this can be expressed as a function of pressure and temperature [6]: where δw is a function of P and Tf, and the value of δw does not exceed 6% of 1/ρw. The temperature of water also affects the kinematic viscosity of water (μ). An empirical formula for kinematic viscosity is [6]: Therefore, the function of density and viscosity of fluid can relate the fluid flow model and heat transport model, and a procedure for coupled simulation of thermal and fluid flow (NPCTF) models in a two-dimension fracture profile was proposed. A detailed flowchart for the proposed simulation method is seen in Figure 3d. When fracture profile information (x,y) is obtained, the fluid flow model is used to evaluate flow rate Q and pressure P, and the heat transport model then predicts the temperature distribution Tf(x). An iterative process is used to solve the fluid flow model and heat transport model by modifying the fluid's physical properties (fluid density ρw and kinematic viscosity μ) using new flowrate, pressure, and temperature before moving to the next step. The iterative process stops once the temperature and flowrate stabilize to a specified tolerance. The tolerances of flowrate and temperature are defined as ||Qj−1 − Qj|| ≤ 0.00001 and ||Tfj−1-Tfj|| ≤ 0.00001, respectively. j is each iteration step in the calculation process. This calculation method would increase calculated efficiency, compared to the previous finite element method or finite difference method.  Barton and Bandis (1985) [41] analyzed the roughness of 136 natural fracture surfaces. According to the roughness of the fracture surface, it is divided into 10 levels, of which values are in the range of 0-20. The roughness index is called the joint roughness coefficient (JRC). The JRC parameter is widely used in rock mechanics and engineering. In order to valid the NPCTF, 5 JRC profiles were selected to calculate temperature and pressure distribution, and those results were compared to numerical results with COMSOL. The JRC values of the selected fracture profiles are 1, 5, 10, 15, and 20, respectively. Table 1 lists the thermal and flow parameters for predicting the temperature in JRC profiles.  Figure 4 shows the comparison between the predicted results of the NPCTF and the calculated results of COMSOL. As illustrated in Figure 4, the calculated results of the NPCTF are consistent with the COMSOL results. With the increase in JRC value, the error becomes larger and a more obvious difference is observed at the nonlinear part of the temperature profile. However, the maximum error value is 12.3%. The reason is that the effect of thermal diffusion is not considered. With the increase in JRC, the thermal diffusion effect is more obvious, and enhances the differences in the nonlinear part of the temperature profile between the proposed NPCTF and COMSOL results. Figure 5 shows the temperature counter in the fracture profile with different JRC. In addition, the results of flow field calculated by proposed NPCTF and COMSOL are compared. The flow field is expressed in terms of transmissivity Ta, which is defined as follows:

Validation of NPCTF
(23) Figure 6 shows the comparison of transmissivity Ta of fractures under different JRC values. The transmissivities predicted by the proposed NPCTF are consistent with the COMSOL results. As the JRC value increases, the error increases, and the maximum error is 1.2%.
Therefore, the proposed NPCTF can predict temperature and flow velocity distribution in rough-walled fractures.

Modelling Setup
To obtain the rock fractures, an intact granite block was split using the Brazilian indirect tensile test into two halves with dimensions of 150 × 150 × 75 mm. This method was also used by Xiong et al. [24] to study fluid flow behavior in rock fractures.
To determine the surface information of the fracture, an advanced 3D CaMega stereotopometric scanning system (BoWeiHengXin Inc., Beijing, China) was used to perform a non-contact 3D scan of the rough fracture surfaces. Then, the coordinates of 22,801 points of the upper half of the fracture were obtained, as well as the lower half fracture. Tatone  [42] introduced a novel method to measure fracture geometry based on fracture surface data. Based on this method, a 3D fracture model was obtained, which is shown in Figure 7a.
Five evenly spaced surface profiles parallel to the flow direction with lengths of 150 mm of the upper and lower half of the fracture were taken into consideration, namely slices X = 15, 45, 75, 105, and 135 mm, as shown in Figure 7a. The profile curves of the five slices were obtained by importing the coordinates (x, z) of every point of each slice to Auto CAD. Correspondingly, different values of JRC are calculated according to Tse and Cruden [43], who proposed a relationship of JRC and the root mean square of the first derivative of profile Z2: Table 2 lists the JRC and mean aperture of the fractures. The heat flow simulations in this work were conducted by solving a 3D Navier-Stokes equation and an energy equation with an advanced COMSOL Multiphysics code. The mesh scheme in the COMSOL mesh module was optimized to reduce the mesh-related effect on flow simulations. The fracture computational mesh was further refined, with denser mesh in the fracture part to capture the details of the flow. The final fracture mesh is shown in Figure 7b. The mean mesh sizes of the rock matrix and fracture are 35.0 and 5.0 μm, respectively. In general, over 3.2 million tetrahedron mesh elements were constructed for the fracture model. The temperature boundary conductions were the upper and lower surface of the rock matrix. No-flux and non-slip boundary conditions were set for the fracture walls. The inflow boundary was assigned a constant mass flow rate and a constant temperature, and the outflow boundary was set as a constant pressure condition.

Temperature and Pressure Distribution Along Fracture
To demonstrate the general thermal transport within the fractures, the temperature counter and distribution are plotted in Figure 8 along five fracture profiles tested in this study using COMSOL simulation results. The calculated parameters are listed in Table 1.
The corresponding injection pressure was 0.01 Pa. As illustrated in Figure 8, the water has a strong trend to attain equilibrium temperature (363.15 K) close to the fracture inlet since the coupling between the water and the rock is strong and, therefore, the heat transport is more intense. Far away from the fracture inlet, the water temperature slowly attains equilibrium temperature. If the fracture is long enough, the temperature at the fracture outlet position would be same as the equilibrium temperature. Additionally, each temperature along the rough-walled fracture is not a smooth curve; the rougher fracture surface leads to more fluctuation for the temperature curve, indicating the fracture roughness affects the temperature distribution, as shown in Figure 8. It can also be seen that as fracture aperture increases, the temperature at which water is absorbed decreases. The difference between the temperature related to the aperture parameter is readily apparent when the change in temperature that occurred between fractures with bm = 1.2 mm and bm = 1.8 mm (the smallest and largest fracture) aperture is examined. A decrease in temperature of nearly 25% was calculated between the flow in these two fractures. This relatively large change in aperture relates well to the large simulated change in temperature. Hence, both fracture roughness and aperture have significant influences on temperature distribution.
The pressure distribution of the fracture is also plotted in Figure 9. As illustrated in Figure 9, from the inlet to the outlet, the fluid pressure decreases from 1 to 0 Pa. Both fracture roughness and aperture have small influences on pressure distribution.

Comparison with Previous Models from the Literature
To demonstrate the robustness of the proposed NPCTF (the calculation process shown in Figure 3) in describing this thermal transport in realistic 3D rock fractures, the deviation D is defined as the relative error in estimating the temperature distribution: 100% (25) where Ti-model is the fracture's local temperature at i node obtained from the NPCTF and Tireal is the fracture local temperature at i node obtained from numerical simulations. The temperature deviations D calculated by Equation (25) for the five fracture profiles are plotted in Figure 10 at different fracture positions. In general, the two temperatures predicted by NPCTF and COMSOL simulations are in close agreement for all fracture profiles, with the deviation D ranging from 0.44% to 8.10% with a mean of 1.51%. A highest mean deviation <D> of 2.57% is observed for the fracture profile with JRC = 16.34. The fracture profile with JRC = 11.96 has the lowest <D> of 0.75% in this case and the range of D is from 0.48% to 3.98%. The root mean squares (RMS) of deviation D were calculated. It is interesting that the fracture profile with the higher JRC has the higher RMS. The highest RMS of 0.0173 is observed for the fracture profile with JRC = 16.34. The fracture profile with JRC = 11.96 has the lowest RMS of 0.00858. The deviation D increases with increases in JRC, also suggesting that the prediction accuracy of the proposed NPCTF is related to fracture roughness. Further assessments of NPCTF were made by comparing the values of D with those obtained from previously proposed heat transport models including the Gehlin and Hellstrom [2] and Zhao models [34]. A steady-state analytical solution describing the temperature distribution in the vertical fracture was developed by Gehlin and Hellstrom [2], as follows: exp (26) where α is the heat transfer coefficient; the value of 900 W m −2 K −1 was used in the research of Shaik et al. [44]. To is rock temperature and Tin is water temperature in the inlet of the fracture. When the diverging angle or converging angle is 0, the proposed heat transport model can be degenerated as in the Zhao model [34]: For comparison, the heat transport model took the place of the Gehlin and Hellstrom and Zhao heat transport models in the NPCTF, when applying the Gehlin and Hellstrom [2] and Zhao [34] models to calculate the temperature distribution in all five fracture cases. The mean temperature deviations <D> (see in Equation (25)) for all five fracture cases by comparing the NPCTF, the Gehlin and Hellstrom model [2], and the Zhao model [34] with COMSOL simulation are plotted in Figure 11. Overall, the highest deviation is observed for the Gehlin and Hellstrom model [2], with a mean of 23.27% and a maximum of 27.81%. The Zhao model [34] deviates by less than 10% of <D> when JRC is less than 13; however, the deviation rises to 12.45%, 12.75%, and 14.16%, respectively, as JRC increases to 16.34. The poorer performance of the two models is due to not taking fracture roughness and aperture variation into consideration, compared to the proposed NPCTF. Therefore, it can be observed that the D from the NPCTF is less than other models. Hence, the proposed NPCTF outperformed other models in different JRC with the lowest <D> of 0.44%.

Sensitivity Analysis
The accuracy of the proposed NPCTF is mainly affected by the fracture geometry (i.e., roughness and aperture variation) and mesh size (one dimensional element). The  Zhao (2014) fracture roughness has been discussed above. Hence, the two factors (aperture and mesh size) need to be tested in the NPCTF.
The fracture profile with maximum roughness (JRC = 16.34 and bm = 1.5 mm) was selected to test, and aperture bm was fixed at 0.6, 0.9, 1.2, 1.5, 1.8, and 2.1 mm by moving the distance of the lower and upper fracture surface. Figure 12 plots the effect of aperture on mean deviation <D> (see in Equation (25)) from NPCTF). As shown in Figure 13, the mean deviation <D> decreases as the value of bm increases. At JRC = 16.34, as the value of bm changes from 0.6 to 2.1 mm, the mean deviation <D> decreases from 1.80 to 0.62. The deviation of cases for aperture less than 1.5 mm is larger than that of cases for an aperture value of 1.8 to 2.1 mm. It is evident that the effect of bm on <D> is larger for a small aperture fracture and smaller for a larger aperture fracture.
In order to analyze the effect of mesh size, the mesh size is set as 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, and 0.9 mm and 1.0, 1.1, and 1.2 mm for the fracture profile with JRC = 16.34 and bm = 1.5 mm, respectively. Figure 13 shows the effect of mesh size on mean deviation <D> (see in Equation (25)) from the NPCTF. As shown in Figure 13, the mean deviation <D> increases as the mesh size increases from 0.6 to 1.2 mm. The maximum deviation for a mesh size value of 1.2 mm is 1.57%. The results tend to converge when the simulation mesh size is less than 0.6 mm. The differences in deviation from using a mesh size of 0.3 and 0.6 mm are approximately 1.38% and 1.41%, respectively. Therefore, this indicates that the NPCTF is not sensitive to the size of mesh defined in the range of 0.2 to 0.8 mm.

Limitations of NPCTF
The main limitation of this proposed NPCTF from the heat transport model is neglecting the heat diffusion term. This assumption can be violated as the diverging or converging angle increases. This heat diffusion effect can only cause bigger deviation in the distance from the fracture inlet to equilibrium temperature position. Furthermore, the deviation from the proposed heat transport model is reasonable when presented a diverging wedge with an angle that is 63.47°. The feature of the diverging and converging angle is mostly considered and reported for rock fractures [36,37,45] and therefore, the validity of the assumption is warranted in most cases. Additionally, the validity of the used flow perturbation solution in the fracture lies in the assumptions that the variation of aperture along the fracture length direction is small. However, Wang et al. [37] discussed that the assumption can meet the geometry of most natural fractures.

Conclusions
We presented an NPCTF for rough-walled fractures, including a fluid flow model and a heat transport model. The fluid flow model describes the relationship of flowrate and pressure in the wedge-shaped cell, based on the perturbation method. The heat transport properties between the flowing fluid and rock fracture wall are captured by the proposed heat transport model, which assumes an identical temperature between the fluid and fracture wall. Both of them connected to each other in a series of connected wedge-shaped fractures. The proposed NPCTF was validated by a commercial COMSOL Multiphysics code in JRC profiles. The temperature relative deviation between them is less than 12.3% and the maximum transmissivity deviation is 1.2%. To examine the performance of the proposed NPCTF, heat flow simulations in real rock fractures with different relative surface roughness, obtained by a 3D stereotopometric scanning system, were performed. It has been shown that the results from the proposed NPCTF agree well with the COMSOL simulation results. In general, the absolute deviation of the proposed model from the simulation results, in terms of the estimated temperature, is in the range of 0.44% to 8.10% with a mean deviation of 1.51% for all cases of the fracture profiles examined. In addition, the NPCTF is sensitive to fracture aperture variation, but not to the size of the mesh. The effect of aperture on deviation of the NPCTF is larger for a small aperture fracture compared with a large aperture fracture.
The proposed NPCTF has good estimation of the heat transport effect of roughwalled fractures. Therefore, in the next work, we would take the proposed NPCTF to use Data Availability Statement: All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest:
The authors declared that they have no conflicts of interest to this work.