Heat transfer characteristics of water flow through a single rock fracture: The combined effects of shear and surface roughness

An accurate understanding of the heat transfer of water through rock fractures is essential for the extraction and utilization of thermal energy from high‐temperature rock masses. A systematic numerical simulation based on the double‐rough‐walled model has been presented to investigate the shear effect on convective heat transfer in rough rock fractures. On the basis of the modified successive random additions algorithm, four different self‐affine surfaces were generated and utilized to establish the 3D double‐rough‐walled fracture models. The fluid flow and heat transfer were simulated by directly solving the Navier–Stokes equation and energy conservation equation, respectively. The combined effects of shear and surface roughness on the heat transfer were investigated. The results show that the heat within rough‐walled fractures transfers preferentially along the main fluid flow channels, and the areas of fast and slow thermal transmission fit well with the high‐ and low‐flow regions, respectively. As shear advances, the heat transfer coefficient firstly increases, then decreases slightly and finally stabilizes within a certain range, in which stabilization occurs earlier in fracture with a larger joint roughness coefficient. The effect of surface roughness on heat transfer shows an opposite trend during shearing. When the shear displacement is small, the enhancement effect of surface roughness that provides larger heat transfer areas dominates the heat transfer. As shear displacement continues to increase, this enhancement effect will be gradually weakened until the decreasing effect that bumps on the rough‐walled surface hinder the fluid flow dominates the heat transfer.


| INTRODUCTION
2][3] Taking the enhanced geothermal system as an example, its essence is the flow and heat transfer process of low-temperature fluid in high-temperature fractured rock mass.A comprehensive understanding of heat transfer of fluid flow through fractured rock is crucial for geothermal reservoir management in terms of optimizing heat extraction and improving the heatrecovery factor of geothermal reservoirs.The process of heated flow through fracture networks is complex owing to the anisotropic geometries of fractures. 4Although the topic has been widely studied in recent years, it is still not completely understood even in the case of a single rough fracture.Flow and thermal modeling for single rough rock fractures is a fundamental problem for assessing the heat transfer through complex fracture networks. 50][11][12][13][14][15][16][17] To describe the fracture morphology, different rough-walled fracture models have been proposed, such as the "single-roughwalled model" in which the lower wall of the fracture is smooth and the upper wall is rough, 18 and the "doublerough-walled model" 19 in which both the upper and the lower wall are rough.Those models are further characterized as two-dimensional (2D) and theredimensional (3D) objects.When modeling the fracture as the 2D plane, the variable apertures are assigned to the fracture elements at a certain interval along the fracture length direction. 6,20,21Since any contact between the walls will hinder the fluid flow, 22 the 2D calculation cannot capture the contact effect, which deviates from the natural state.Wang et al. 23 proposed a 3D doublerough-walled model by building a rough fracture wall and then adding a normally distributed aperture to each element to form another rough wall.The model successfully takes into account the anisotropic aperture distribution, but fails to reflect the void space induced by shear.Another commonly used 3D double-rough-walled model is to create two well-matched rough surfaces, and then artificially move the upper and lower walls of the fracture by assigning shear displacement and normal displacement. 24This model simplifies the correlation between the normal displacement and shear displacement, which makes it difficult to describe shear-induced dilation.
Due to natural events and human activities, rock fractures are usually subjected to in situ stresses containing normal stress and shear stress, which alters the geometric morphology and aperture of fractures, thereby affecting the hydraulic and heat transfer properties of fractures.2][33] During the fracture closure under normal stress, the fracture aperture decreases and the contact area between two surface walls increases, which results in a reduction of fluid flow thus affecting the heat transport [34][35][36][37] built a series of 3D fracture models by measuring the fracture morphologies under different normal stresses, and the effect of fracture geometrical alterations induced by changing stresses on the heat transfer processes was investigated.They proposed that the anisotropic aperture field leads to heterogeneous distributions of streamlines and water temperatures, and such heterogeneity becomes more obvious with increasing fracture closure.Sun et al. 38 proposed that the normal stresses induced aperture closures would result in the delay of thermal breakthrough compared with the unstressed boundary case.According to Huang et al. 39 under a high hydraulic gradient, increasing stress ratio may lead to the emergence of preferential flow paths and a higher network connectivity, which inhibits sufficient heat exchange between fluids and rock matrix, leading to a low heat extraction.The previous studies mentioned above have explored the fluid flow and heat transfer process under normal stress, but the fluid flow and heat transfer characteristics during shearing have not been thoroughly understood.
In light of the abovementioned issues, it is important to clarify the fluid flow and heat transfer characteristics under shear stress in complex high-temperature rock environments.In this study, a series of geologically realistic fracture models were established, in which the rough walls of fracture were generated using a modified successive random additions (SRAs) algorithm, and aperture distributions under the shear process were determined based on a mechanical model.The fluid flow and heat transfer processes of rock fractures were simulated by solving the Navier-Stokes equation and the energy conservation equation, using the numerical analysis software COMSOL Multiphysics.The effects of shear and surface roughness on fluid flow and heat transfer were systematically investigated.The evolutions of temperature of outlet water, energy extraction efficiency, and heat transfer coefficient with shear displacement were quantified.

| Generation of rough fracture surface
The surface of natural rock fractures is usually rough and follows a self-affine fractal distribution, which can be modeled using the fractional Brownian motion (fBm). 40urrently, many methods have been proposed to generate the fracture surfaces using fBm, such as SRA, Weierstrass Mandelbrot function randomization, and Fourier transform. 41,42In this study, the SRA algorithm is used to simulate the heterogeneity of fracture surfaces due to its simplicity and efficient.
In the 2D fBm, the continuous single-valued function Z(x, y) can be used to represent the aperture asperity height of a fracture wall.The stationary increment [Z(x + l x , y + l y ) − Z(x, y)] over the distance (lag) l has a Gaussian distribution with a mean of zero and a variance of σ 2 . 43,44The statistical self-affine property of fBm increment can be expressed as (2) where 〈•〉 represents the mathematical expectation, H is the roughness exponent or Hurst exponent varying from 0 to 1 and related to the 3D fractal dimension (D f ) by 23,45 r is a constant, and σ 2 is the variance defined as (3) Liu et al. 46 have developed a modified SAR algorithm to generate 3D self-affine surfaces of fractures, which solves the problem of generating random fractal distributions with poor correlation in traditional algorithms.In this study, a series of 3D self-affine surfaces with H = 0.50, 0.55, 0.60, and 0.65, respectively, were generated using the modified SAR algorithm proposed by Liu et al. 46 The detailed steps of the modified SRA algorithm have been displayed in the previous study of Liu et al., 46 and thus we refer the author to them for more details.A larger H results in smoother fracture surfaces, and vice versa.The surfaces are 256 mm × 256 mm in size, having 513 points with an interval of 0.5 mm along the x-and ydirections, respectively (Figure 1).
The joint roughness coefficient (JRC) for the fracture profile is calculated using the following formula proposed by Tse and Cruden 47 : where Z 2 is the root mean square of the first derivative of the profile, and can be expressed in the discrete form where x i and z i represent the coordinates of the fracture surface profile, and N t is the number of sampling points along the length of a fracture.The average of JRC values for a series of profiles along a surface in the direction parallel to the shear direction is calculated 48

| Aperture assignment
The fracture aperture is closely related to permeability and heat transport, so it is very important to accurately simulate the variation of fracture aperture during shearing.In this study, it is assumed that the two walls of the fracture are well mated with an initial aperture of zero before shear begins.When the modified SRA algorithm is used to generate a fracture surface, the fracture aperture b at the shear displacement u s can be expressed by the function Z(x, y): where u v is the normal displacement at u s .Indraratna et al. 49 proposed the following combined equations to calculate u v : where v̇is the dilation rate, u s-peak is the peak shear displacement, c 0 is the ratio of u s /u s-peak at which dilation begins, c 1 and c 2 are two decay constants, and vṗ eak is peak dilation rate that can be calculated according to where S n is the normal boundary stiffness, α, β, and λ are parameters that can be determined using the following equations: where JCS is the compressive strength of the joint surface, k ni is the initial fracture normal stiffness at zero normal stress, V m is the maximum closure of the fracture, σ n0 is the initial normal stress, and M is the damage coefficient that equals to 1 or 2 depending on the magnitude of normal stress. 50 small initial normal stress of 1 MPa is applied, so that the asperity degradation or surface damage is not considered in the present primary study.Many previous studies have investigated the shear effect on the fluid flow through rock fractures using laboratory experiments and numerical methods. 49,51,52Referring to the shear displacement considered in those studies, the present study selects a slightly larger range of 0-20 mm, and picked up several shear displacement intervals of 1, 2, 3, 5, 10, 15, and 20 mm.
Table 1 summarizes the parameters used in the simulation.Since the nominal contact between the upper and lower surfaces of fractures decreases as the shear displacement increases, the fracture model with a crosssectional area of 200 mm × 100 mm in the xy-plane was extracted from the central part of the original model after shear to maintain a consistent analyzed model size.
Figure 2 shows the space distribution of the aperture from JRC4.
The aperture distributions of fracture models with JRC1-JRC4 under different shear displacements were calculated as mentioned above, and the results are displayed in Figure 3.It can be seen that the aperture distribution is well approximated by the normal distribution.Larger shear displacement results in larger variance and mean value, and thereafter, the larger void spaces between the upper and lower fracture surfaces.This indicates that the aperture field becomes more anisotropic and heterogeneous as the shear displacement increases.
The part with an aperture less than 0 is defined as the contact area between the upper and lower surfaces of fracture, which does not participate in fluid flow.At the laboratory scale, when the shear displacement of a wellmated fracture is small, it is difficult to form a connected flow channel from the inlet to the outlet of the fracture.As with many previous studies, 9,37,53 to ensure at least one connected flow path, the upper surface of fracture after each shear displacement is slightly shifted in the vertical direction.In this study, the lifting height is selected to be 0.05 mm, which is the minimum vertical shifting height to ensure the fracture is permeable at all shear displacements.

| NUMERICAL MODELING AND SIMULATION
where u is the velocity vector (m/s), μ and ρ are the viscosity (Pa s) and density (kg/m 3 ) of the fluid, respectively, and P is the fluid pressure (Pa).
The continuity equation can be expressed as

| Governing equation of heat transfer
Since the variation of fluid temperature is small in the simulations, it is assumed that no fluid phase changes during the flow process, and the density and viscosity of fluid remain constant.Those assumptions may give rise to some errors when evaluating the properties of heat transfer, which is neglected in this study.In further study, we will extend our investigations to incorporate the fluid phase transition and viscoelastic effects.
The two commonly used continuum mechanical approaches describing heat transport in media include local thermal equilibrium (LTE) and local thermal nonequilibrium (LTNE).The LTE theory assumes that the temperature of the fracture surface is consistent with the fluid at that location, and the temperature distribution of both rock and fluid can be expressed by one equation to describe the entire system.In the LTNE | 565 formulation, each phase has a distinct temperature at material points and two temperature fields must be considered.Compared with the LTE, the LTNE is more precise to describe the heat transfer between the fracture surface and the fluid; whereas, it needs the explicit parameterization of the heat exchange term as well as a large computational demand.5][56] Heat transport in a rock matrix is controlled by conduction, and the energy conservation equation of rock under steady conditions is expressed as where T r and K r are the temperature (K) and the thermal conductivity of the rock (W/m K), respectively.Heat transport in fluid is controlled by convection and conduction.The energy conservation equation of the fluid in the fracture under steady conditions is expressed as where C p is the specific heat capacity of the fluid (J/kg K), and T f and K f are the temperature (K) and the thermal conductivity of the fluid (W/m K), respectively.

| Numerical simulation setup
The finite element software COMSOL Multiphysics was used to solve the Navier-Stokes equation and the energy conservation equation to simulate the fluid flow and heat transport process in a single fracture.On the basis of the 3D self-affine surface data generated in Section 2.2, the high-precision parameterized surface function in COM-SOL was used to reproduce the real 3D fracture geometries and irregular spatial structure sandwiched between the two surfaces.The 3D numerical model is 200 mm × 100 mm × 100 mm in size as shown in Figure 4.
When solving the equations for fluid flow and heat transport, many previous studies [57][58][59] usually assign very small aperture values to the contact area to avoid solving an ill-formed matrix.As a result, there are still some fluid flows inside the contact areas, which are physically nonrealistic.In this study, the contact regions are regarded as zero aperture elements that do not participate in the flow and heat transfer of the fluid.This will result in more precise flow fields and enable more accurate characterization of the effects of fracture morphology and shear on heat transfer.The tetrahedral elements are utilized to mesh the 3D geometric model, which allows for an accurate discretization around the contact area as shown in Figure 5.
Fluid flow is assumed to occur in fractures with impermeable rock matrix.As shown in Figure 4, a constant flow rate of U in = 0.02 m/s is applied on the inlet boundary, and a zero pressure of P = 0 is applied on the outlet boundary.Note that when changing the flow rate, more complicated flow patterns may occur and thereby affect the heat transfer process, which is beyond the scope of the present study and will be investigated in further studies.The upper and lower boundaries of the fracture are set as no-flow and no-slip walls.For the boundary conditions of heat transport modeling, the initial temperature of both rock and water is set as 373.15K.The cold water of 293.15K is injected along the inflow boundary.The remaining boundaries of the entire model are set as adiabatic boundaries.The total duration of this simulation is 900 s with a time step of 10 s, after which the transient temperature of the fracture and the fluid flow rate remain stable.

| Model validation
This study will use the experimental results obtained by Zhao and Tso 60 as the basis for model validation.More details about the experiment can be found in Zhao and Tso. 60Referring to experiments, a 2D numerical model as shown in Figure 6 was established.The entire model is divided into three layers along the y-direction, in which the upper and lower layers are impermeable rocks, and the middle part is flowing fluid in the rock fractures.The model has a length of 102 mm in the x-direction and a width of 51 mm in the y-direction, which is consistent with the sample size in the study of Zhao and Tso. 60The fluid moves along the x-direction, and the y-direction represents the thickness of the rock.The parameters of rock and fluid used in numerical simulation are shown in Table 2.
The established 2D model was calculated referring to the flow rate and injection temperature used in the experiments of Zhao and Tso 60 as shown in Table 3, where T ext is the temperature of the external surface of the rock sample, and T in and T out represent the inlet and outlet water temperatures, respectively.The simulated outlet water temperatures were compared with the experimental results, and the goodness between the two data sets was quantified by the relative error D, which is defined as Where T s and T e represent the numerical and experimental results, respectively.From Table 3 it can be seen that the simulated temperature is slightly higher than the experimental result, but the relative error is restricted within 3.37%, which verifies the validity of the numerical simulation method.
F I G U R E 4 Initial and boundary conditions for a numerical model of rough fractures.

| Spatial flow and temperature distribution
Fluid flow and heat transfer in rock fractures with different JRCs are simulated, and without loss of generality, Figure 7 only illustrates the results of the fracture model with JRC2 = 7.49 where the influence of shear displacement on the aperture distribution, temperature field, and velocity field can be observed.The white regions shown in Figure 7 represent the contact areas.As shown in Figure 7A,B, a few flow channels that follow a tortuous path bypassing plenty of small contact obstacles appear in the void areas, which is known as the channeling effect. 61As the shear moves on, an increase in aperture and a decrease in contact areas result in more continuous and obvious flow paths.Following those main flow channels, the temperature near the inlet boundary decreases with time and moves heterogeneously towards the outlet direction, as shown in Figure 7C.Such heterogeneity of heat distribution becomes more obvious with increasing the shear displacement.Since the inlet velocity is constant, the increase in the aperture results in a larger volumetric T out (°C) flow rate, which increases the total amount of heat absorbed by the water flowing in the fracture and decreases the heat stored in the surrounding rock.The temperature distribution is sufficiently relevant to the surface morphology of fractures, so the velocity and temperature fields in fractures with different JRCs are discrepant.As shown in Figure 8, the flow path and temperature distribution have similar nonuniformity characteristics, which fluctuate with fracture morphology.In the high-flow rate regions, the heat transport process is dominated by the heat convection, resulting in an obvious decrement of temperature along flow streamlines.In contrast, the heat transfer process in flow-stagnated regions has a small temperature change.Under the same shear displacement, a more sufficient heat transfer between the fluid and the rock is observed in the fracture with a larger JRC, accompanied by a faster temperature drop and lower outlet temperature.This is because the roughness of the fracture surfaces increases the length of the flow path, which contributes to the expansion of flow areas and heat transfer efficiency.

| Temperature of outlet water
Figure 9 shows the change of outlet temperature T out with time for fractures with different JRCs.The results indicate that for all cases T out decreases nonlinearly with time, in which a more obvious decrement is observed at a larger shear displacement.Once shearing occurs, even a small shear displacement will cause a significant decrement in T out .However, when the shear displacement increases to a certain value, the enhancement of convective heat transfer induced by shear is no longer obvious.Therefore, a certain degree of hydraulic shear can effectively promote heat transfer.
To quantify the temperature change during the heat transfer process, ΔT that is defined as the difference F I G U R E 7 In JRC2, the diagrams of (A) the space distribution of apertures, (B) principal velocity field, and (C) fluid temperature fields at t = 900 s when the shear displacements are 5, 10, and 20 mm, respectively.JRC, joint roughness coefficient.between the T out at 0 and 900 s is calculated.Figure 10 shows that as shear displacement increases, ΔT increases rapidly first and then gradually becomes stable after the shear displacement of 10 mm.Surface roughness also plays a significant role in the change of ΔT, in which a larger JRC results in a higher ΔT.This indicates that the heat exchange between water and rock will be more adequate in fractures with rougher surfaces.On the other hand, the roughness will increase the bumps and recesses in fracture surface, which would hinder the fluid flow.As a result, ΔT slightly decreases as JRC increases from JRC3 to JRC4, indicating a restricted effect of roughness on promoting the heat transfer.

| Energy extraction efficiency
To account for the energy extraction efficiency, the instantaneous heat extraction Q I defined as follows is calculated: where c p is the specific heat capacity of water (J/kg K), u out is the outlet flow rate (m/s), and b is the average fracture aperture.Figure 11 indicates that when the heat transfer process occurs, the Q I increases rapidly to a peak value.A larger shear displacement always implies a faster increasing rate of Q I and a higher peak value.The peak Q I of JRC1-JRC4 when u s = 20 mm is 1.18, 7.64, 38.96, and 30.26J, respectively.It can be seen that when the surface roughness increases from JRC1 to JRC3, the peak Q I increases, while when the JRC continues to increase to JRC4, the peak Q I decreases slightly, which indicates a restricted effect of roughness on increasing the Q I .As the heat transfer process continues, the Q I decreases due to the reduction in outlet temperature caused by the heat exchange and eventually stabilizes.
The cumulative heat extraction Q c over the simulation time t is calculated as follows: F I G U R E 8 Diagrams of (A) the space distribution of aperture, (B) principal velocity field, and (C) fluid temperature fields from JRC1 to JRC4 for a shear displacement of 10 mm (t = 900 s).JRC, joint roughness coefficient.
Figure 12 presents that the evolution of Q c with t for fractures of different JRCs and u s shows a similar increasing tendency, in which a higher fracture roughness implies a greater increasing rate of Q c .Specifically, when the shear displacement is 20 mm, the Q c of fractures with JRC1-JRC4 at 900 s are 0.244 × 10 5 , 1.09 × 10 5 , 2.62 × 10 5 , and 2.09 × 10 5 J, respectively.This indicates a restricted effect of roughness on increasing Q c , which is consistent with the effect on the Q I shown in Figure 11.The curve increases faster with the increasing shear displacement, due to the combined effect of shearenhanced outlet volumetric flow and the shear-reduced outlet temperature.When the shear displacement is 20 mm, the Q c of JRC1-JRC4 compared with that of 2 mm shear displacement increases by 509%, 2533%, 764%, and 1161%, respectively.It can be seen that the increase in shear displacement can indeed greatly improve the heat-recovery efficiency of the thermal reservoir system.
The results of this study may provide a basis for fieldscale modeling of thermal reservoirs.The increasing flow flux and heat extraction efficiency induced by shear should be fully recognized and utilized.At the same time, the acceleration of heat extraction will also shorten the service life of a thermal reservoir system.Therefore, balancing the shear-enhanced heat extraction efficiency and heat reservoir lifespan would be significant to the sustainable development and utilization of geothermal energy.

| Heat transfer coefficient
The heat transfer coefficient can be used as an important parameter to evaluate the heat transfer ability of a single fractured rock mass.In this study, the heat transfer coefficient is calculated using the following equation proposed by Bai et al., 62 which assumes that the temperature along the flow direction is a linear function: where h is the heat transfer coefficient (W/(m 2 K)), c p,w is the specific heat capacity of water at constant pressure (J/(kg K)), ρ w is the density of water (kg/m 3 ), q v is the volumetric flow rate (m 3 /s), d is the diameter of the rock (m), L is the length of the rock (m), and T c is the temperature of the external surface of the rock (K).
As shown in Figure 13, the h of four fractures all increases rapidly before u s = 3 mm, indicating that the h will increase substantially once shear occurs.This implies that a small shear displacement can greatly promote the convective heat transfer.As the shear F I G U R E 10 Comparison of ΔT under different JRCs and u s .ΔT that is defined as the difference between the T out at 0 and 900 s.JRC, joint roughness coefficient.continues, the h decreases slightly and finally stabilizes within a certain range, in which the h stabilizes earlier in fracture with a larger JRC.The coupling effect of shearenhanced volumetric flow rate and shear-reduced outlet temperature controls the heat transport.When shear occurs, the rapidly increased aperture leads to a higher volumetric flow rate and a greater total amount of water that absorbs heat, thus a markedly increased h.However, as the fracture aperture increases with continuously increasing shearing, both the outlet temperature of rock fracture and the heat energy stored in the rock decreases to some extent, resulting in a small reduction in h.
Finally, when the influence between volumetric flow rate and outlet temperature reaches equilibrium, the h stabilizes within a certain range.Figure 13 indicates that the effect of surface roughness on h shows an opposite trend as the shear displacement increases, in which larger JRC results in a larger h when u s < 3 mm but a smaller h when u s > 15 mm.This is induced by the two competing effects of the surface roughness on the heat transfer during shear: increasing surface roughness can lead to larger heat transfer areas which would promote the heat transfer; besides, the tortuosity will increase due to the bumps and recesses on the fracture surface with increasing roughness, which would hinder the fluid flow and slow down heat transport.When u s < 3 mm, the small aperture restricts the amount of fluid passing through the fracture, so that the enhancement effect of surface roughness dominates the heat transport.This enhancement effect will be gradually weakened as the aperture increases with increasing the shear displacement.When u s > 15 mm, the reduction effect of surface roughness dominates the heat transport, thus resulting in a smaller h for a larger JRC.In general, surface roughness has a promoting effect on heat transfer when the shear displacement is small.With the shear displacement increasing, this enhancement effect will be gradually weakened.As the shear displacement is large enough, the surface roughness effect becomes less prominent.It should be noted that the h of JRC3 is smaller than that of JRC4 for a larger u s .This is because when the fracture hydraulic aperture is above a certain threshold, the effect of surface roughness on heat transport decreases with the increase of volumetric flow. 37The larger volumetric flow rate of JRC4 increases the h and offsets the decrement of h induced by surface roughness, thereby resulting in a larger h than JRC3.The inflection point of the heat transfer coefficient at u s = 3 and 15 mm in Figure 13 only shows the trend of the competing outcome of the two effects during shearing.Further calculations are needed to determine the exact critical value, which may depend on surface roughness, shear displacement, as well as inlet velocity.In the future, when a broader scope of fracture morphology, shear displacement, and inlet velocity are considered, we will try to give a universal quantification of this critical value.

| CONCLUSION
In this study, a modified SRA algorithm was used to generate a series of geologically realistic fracture rough walls.The aperture distributions in those rough-walled fractures under the shear process were determined based on a mechanical model.The fluid flow and heat transport processes of rock fractures were simulated by solving the Navier-Stokes equation and energy conservation equation using the COMSOL Multiphysics.The effects of shear and surface roughness on the heat transfer characteristics of water flow through rough fractures were investigated.The following conclusions can be obtained: (1) The heat within the rough-walled fractures transports preferentially along the main fluid flow channels, and the areas of fast-and slowtransmission fit well with the high-and low-flow regions, respectively.(2) The heat exchange between water and rock will be more adequate in fractures with rougher surfaces.As the continuous increase of JRC, the phenomenon of promoting heat transfer is no longer obvious, indicating a restricted effect of surface roughness on the heat transfer intensity.This is the result of the combined effect of rough surfaces increasing heat transfer area and bumps hindering effective heat transport.(3) When the surface roughness remains constant, the heat transfer coefficient firstly increases with shear, then decreases slightly, and finally stabilizes within a certain range, in which the heat transfer coefficient stabilizes earlier in the fracture with a larger JRC.As the shear enhances the hydraulic conductivity of rock fractures, the instantaneous heat extraction and cumulative heat extraction increase with the increase of shear displacement.(4) The effect of surface roughness on heat transfer shows an opposite trend as the shear displacement increases.When the shear displacement is small, the enhancement effect of surface roughness that increases the flow path and provides larger heat transfer areas dominates the heat transfer.With the increase of shear displacement, this enhancement effect will be gradually weakened until the decreasing effect that bumps and recesses on the roughwalled fracture surface hinder the fluid flow dominates the heat transport.
Flow and thermal modeling for single rough rock fractures is a fundamental problem for assessing heat transfer through complex fracture networks.In the future, we will extend our study to complex fracture networks, and quantify how the shear and surface roughness impact the heat transfer characteristics, which are potentially correlated to the geometrical properties of fracture networks (e.g., fracture density, orientation, length, etc.).

3. 1 |
Governing equations 3.1.1| Governing equation of fluid flow For stable flow, the motion of the fluid through a single fracture is controlled by the Navier-Stokes equation, which can be expressed as

3
Aperture distributions of fractures: (A)-(D) correspond to the line chart of the aperture from JRC1 to JRC4.b ̅ is the average aperture; σ b is the variance.JRC, joint roughness coefficient.

F 3
I G U R E 5 Tetrahedral mesh of the FEM model.The mesh of two rock blocks is not shown.FEM, finite element method.F I G U R E 6 A two-dimensional concept model of heated flow through in a rock fracture.T A B L E 2 Input parameters for numerical simulations.Parameter K (W/m K) ρ (kg/m 3 ) C p (J/kg K) Experimental results and numerical simulation predictions for heated flow through in a rock fracture.

F
I G U R E 11 Relationship between Q I and u s for different JRCs.(A) JRC1, (B) JRC2, (C) JRC3, and (D) JRC4.JRC, joint roughness coefficient.

F I G U R E 12
Relationship between Q c and u s for different JRCs.(A) JRC1, (B) JRC2, (C) JRC3, and (D) JRC4.JRC, joint roughness coefficient.
Geometrical properties of the fracture model and input parameters for aperture assignment.
T A B L E 1 F I G U R E 2 An example of the void geometry of fracture model with JRC4.The color intensity indicates the magnitude of asperity height.JRC, joint roughness coefficient.WANG ET AL.