A fully coupled seepage–heat transfer model including a dynamic heat transfer coefficient in fractured rock sample with a single fissure

Abstract Conventional seepage–heat transfer models in simulating the heat transfer between fluid and rock in fractures mainly involve one-way coupling and do not consider the influence of temperature on the seepage. Moreover, it is an enormous challenge to define parameters of the explicit heat transfer between the rock and fluid. In order to resolve these shortcomings, a two-way fully coupled model of the seepage–heat transfer in the fractured rock was established in the present study. Based on the original geometric structure of the experimental device, combined with the actual engineering scale, the local dynamic heat transfer coefficient of the fractured rock was established, which is related to the fracture aperture, flow velocity and thermal parameters. Then, the proposed model was verified through the experiment and excellent agreement was achieved in this regard. It was found that the dynamic heat transfer coefficient changes the temperature distribution of fluid and rock in the original static heat transfer coefficient fracture system. The proposed model simplifies the parameters required for the calculation of the heat transfer coefficient. These parameters are related to characteristic variables, such as velocity and rock temperature, and can be simply obtained from standard laboratory tests.


Introduction
With the continuous increase in energy demand in the past few decades, large-scale rock mass projects such as oil and gas extractions, nuclear waste disposal, underground oil and gas storage, deep geothermal exploration, and tunnel excavation have emerged worldwide. These projects usually involve seepage and heat transfer in the fractured rock masses. Studies show that the influence of seepage on the temperature field of the fractured rock mass mainly originates from the fluid movement, which promotes the heat exchange between the rock matrix and the fluid. In the convective heat transfer mechanism, the fluid operates as the heat carrier, so that the inflow and outflow of the fluid affect the temperature distribution in the rock mass (Bataille et al. 2006;Jiang et al. 2014;Chen et al. 2018;Bauer et al. 2019;Ndikubwimana et al. 2020;Sun et al. 2020;Xu et al. 2020).
Reviewing the literature indicates that numerous experimental and numerical investigations have been carried out on the flows and heat transfer in fractured rocks. However, most of them are about the seepage in rock fractures, while only a few experiments have been conducted on the coupled seepage and heat transfer phenomenon. Experimental methods in this field can be mainly divided into two categories: laboratory experiments and field experiments. Considering the massive scale of fractured rock, it is often a challenge to apply the results of laboratory experiments to practical projects. Accordingly, laboratory investigations have mainly focused on a single fissure, and the fracture is usually prepared by splitting cylindrical granite in the Brazilian split test (Zhao and Tso 1993;Lu and Xiang 2012;Bai et al. 2017;Li et al. 2017). In other words, there are few experimental studies on the coupling of seepage and heat transfer in a single fissure (Rashad et al. 2014;Wang et al. 2017). Zhao conducted hydraulic-thermal experiments on granite samples to study the influence of the temperature field on the permeability of fractures and the convective heat transfer coefficient between the fracture water and the rock matrix. Then the obtained results were compared with empirical equations based on the heat transfer theory, and it was found that the convective heat transfer coefficient is related to the seepage velocity of the fractured water (Zhao 1994). The conclusion that the convective heat transfer coefficient was related to the seepage velocity of fractured water was compared with empirical equations on the theory of heat transfer. Field experiments have been mainly based on hydraulic experiments and tracer experiments. Many researchers have also performed numerical simulations on the seepage and heat transfer process of fractured rocks, and studied the thermal-hydraulic coupling mechanism and parameter sensitivity (Cui et al. 2016;Jiang et al. 2018;Bauer et al. 2019;Han et al. 2020). Some numerical studies have shown that the heat transfer during the seepage and heat transfer process of the rock fracture system plays an important role (Ewing and Wang 2001;Abbasi et al. 2017).
Many investigations have been carried out on the seepage-heat transfer in porous media and variations of coupling parameters (Junrui 2002). However, the coupling of seepage and temperature fields (hydrothermal coupling) in fractures has been addressed in only a few studies (Gringarten et al. 1975). Aiming at the problem of fracture water flow and heat transfer in engineering applications, a simplified mathematical model was established to study the water flow and heat transfer in multiple groups of parallel fractures. In this model, the geothermal gradient, was considered and an analytical solution was obtained for the temperature field of the fracture water with only a single fissure in the rock mass (Yao et al. 2018;Xu et al. 2020b). Based on the local thermal nonequilibrium theory, a mathematical model consisting of the thermal-hydrological-mechanical (THM) coupling process and an ideal 3D-EGS numerical model was established to simulate the heat generation of EGS and the distributions of pressure, temperature, stress and deformation in a geothermal reservoir (Sun et al. 2020). Then a mathematical model of the thermal-hydrological (TH) process in fractured rock masses was developed based on the thermal nonequilibrium theory. The accuracy of the established model was then verified through analytical solutions. It is worth noting that most established models are one-way couplings of seepage and heat transfer that do not consider the influence of temperature on the seepage field. Investigating the seepage-stress coupling effect of the project shows that the seepage effect in the fissure has a fatal effect on the stability of the surrounding rock (Junrui and Yanqing 2001;Junrui et al. 2004Junrui et al. , 2020. It should be noted that the overall form of the fractured rock seepage-heat transfer model has been obtained and the main challenge is to study the coupling parameters. The heat transfer coefficient can be used to analyze and evaluate the heat transfer capacity between the fluid and the fracture surface. Meanwhile, it is an important parameter of the fracture seepage-heat transfer coupling model. Currently, different methods have been proposed to estimate the heat transfer coefficient. Studies show that different calculation methods to obtain the seepage-heat transfer model can affect may lead to different estimates of the outlet temperature (Bai, Zheng and Nakayama 2019). The heat transfer coefficient in rock fractures is mainly based on two assumptions. The first assumption is that the temperature is the same between the fluid and the fracture surface (Zhao and Tso 1993;Zhang et al. 2015). The second assumption is that the heat transfer coefficient is constant along the fracture surface so that an average value is considered in the calculations (Shaik et al. 2011). However, actual experimental data show that, in practical applications, the heat transfer coefficient in the fracture is correlated with other variables and parameters in the model such as flow velocity, fracture aperture, and temperature. Therefore, the latter assumption is rarely used in geothermal simulations.
In the present study, a locally defined dynamic heat transfer coefficient equation was proposed. Then a two-way fully coupled equation of seepage-heat transfer in fractured rocks was established to resolve the shortcomings of conventional models. The established model is expected to achieve accurate results in simulating the heat transfer of the rock fracture system in the laboratory and meet the evolving heat transfer scenarios in practical and engineering scales. A numerical simulation was carried out to evaluate the accuracy of the proposed model using the experimental data (Zhao and Tso 1993). The parameters and heat transfer efficiency involved in the dynamic heat transfer coefficient are analyzed to prove the scalability of the model.

Basic assumptions
Based on certain facts, the following assumptions are made in this paper.
(1) Seepage only occurs in the fracture and conforms to the cubic law.
(2) The fracture morphology is ignored in the research.
(3) The permeability of the rock and the chemical reaction with the fluid are ignored. (4) The heat conduction in the rock matrix and fracture satisfies Fourier's law. The heat exchange between the rock matrix and the fracture satisfies Newton's law of cooling.

The governing equations for seepage field
When the flow velocity is small, it can be simulated with Darcy's law. Based on Darcy's linear flow, the fractured rock can be simplified to a flat plate model. Moreover, the movement of fluid in the fracture follows N-S equations. Governing equations for the steady-state flow can be expressed as the following: where t is time, q is the density, u is the flow velocity, F is volume force, P is the pressure, and l is the dynamic viscosity, r is Laplacian.
where the first and second terms on the left side of Eq. (1) reflect the velocity acceleration, and the convective acceleration respectively. Moreover, the three terms on the right side are volumetric force, pressure gradient and viscosity, respectively.
For stable laminar flow, @u=@t ¼ 0, and Eq. (1) becomes Usually, temperature and dynamic viscosity are considered constants with fixed values. Therefore, in the seepage-heat transfer coupling model, the seepage velocity field is used to unilaterally affect the temperature field, while the effect of the temperature field on the seepage field is neglected (Shaik et al. 2011;Cao et al. 2019;Ma et al. 2020). It is worth noting that the influence of the temperature on the seepage field can be indirectly reflected by the variations of the fluid properties. To realize the full hydrothermal coupling, the seepage field equation can be expressed as follows: where fluid density and dynamic viscosity are functions of temperature.

The governing equations for temperature field
The governing equation of the temperature field is derived from the two parts of the complete rock and the fracture.
(1) Rock matrix temperature field equation The rock is dense and impermeable, so the heat transfer occurring in the rock matrix is mainly heat conduction. The energy conservation equation in steady state rock is expressed as q r c r @T r @t þ q f c f uDT r þ rq ¼ 0 (4) According to Fourier heat transfer q ¼ Àk r DT r (5) Therefore, the final rock matrix heat transfer equation is where q r is the density of the rock, c r is the specific heat capacity of the rock, t is the time, and T r is the rock matrix temperature, k r is thermal conductivity of the rock, q is specific heat flux.
(2) Rock fracture temperature field equation There are three forms of heat transfer in fluid flow in the fracture. The first form is the convective heat transfer between the fluid and fluid caused by the relative motion of water molecules. The convective heat transfer flow is Q 1 where q f is the density of fluid, c f is the specific heat capacity of the fluid, and T f is the fluid temperature. The second form is heat conduction between fluids, which causes heat transfer flow Q 2 where k f is the thermal conductivity of the fluid. The third form is the convective heat transfer between the fluid and rock fracture surface. According to Newton's cooling law, the convection heat transfer coefficient is used to calculate the heat transfer between the fluid and rock fracture surface. The heat transfer flow is Q 3 where h is the convection heat transfer coefficient and A is the heat exchange area. According to the conservation of energy, the heat required for the temperature change in the rock fracture is equal to the algebraic sum of the heat in and out; therefore,

Physical and thermodynamic parameters
(1) Basic physical parameters The viscosity coefficient of liquid water is a function of temperature.
l ¼ 2:1 Â 10 À6 exp 1808:5 273:5 þ T The density of water changes with temperature and can be expressed as a function of temperature.
(2) Effective heat transfer coefficient The heat transfer coefficient is an important parameter of the convective heat exchange mechanism between fluid and fracture surfaces. Equation (9) indicates that the heat exchange can be calculated using the transfer term based on the temperature difference. Zhao and Tso (1993) and He et al. (2019) proposed the following empirical equations to obtain the heat transfer coefficient. h h where b is width of the fracture, T 0 is the constant temperature at the outer wall surface of specimen, T out and T in are the temperatures of the inlet water and outlet water, L e is the equivalent thickness of half the fracture block, L is the length of the specimen, d is the diameter of the specimen. It should be indicated that Eq. 14 can be derived from one-dimensional heat conduction, while Eq. 15 is based on the heat transfer coefficient in the two-dimensional heat conduction. Eqs. 14 and 15 are both the average heat transfer coefficients obtained from the inlet and outlet temperatures. It is assumed that the temperature along the flow direction has a linear variation and high-order variations are neglected. In the experiment, it is an enormous challenge to measure the temperature of the flow along the fracture direction. Theoretical analyses reveal that the water temperature in the fracture is nonlinearly distributed. Accordingly, the water temperature in the fracture will be underestimated in the linear assumption (Bai et al. 2017;Dirker et al. 2017).

(a) Dynamic heat transfer coefficient based on the boundary layer theory (Model I)
Based on the boundary layer theory, the heat transfer coefficient can be expressed as follows (Lienhard and Catton 1986).
where Re is the Nusselt number, Pr is the Prandtl number, and the expression is as follows: Pr The density and dynamic viscosity of water are functions of the temperature. Accordingly, variations of the fluid temperature affect the proposed heat transfer coefficient which is based on the boundary layer theory.
Substituting Eqs. (12) and (13) into Eq. (19) reveal that the heat transfer coefficient h varies with temperature. This phenomenon is also shown in Figure 1. Figure 2 indicates that the heat transfer coefficient is not constant and changes with temperature. The heat transfer temperature of the flowing water and the inner surface of the fracture change and the process is accompanied by variations of the heat transfer coefficient. This method is used as a fully coupled seepage-heat transfer coupling model I.

(b) Dynamic heat transfer coefficient based on local heat transfer (Model II)
Zhao assumed a constant heat transfer coefficient along the fractures and proposed the following expression for a uniform heat transfer coefficient based on nonlinear changes in water temperature under steady-state condition (Zhao and Geotechnics 2014): Equation (20) shows that the water temperature is known along the fracture surface as T f ðx ¼ LÞ ¼ T out : However, this equation is based on the average heat transfer coefficient calculated by using the difference between the outlet and inlet water temperatures and cannot obtain the local heat transfer coefficient. T 0 denotes the dynamic heat transfer coefficient along the fracture surface, so Eqs. (20) (14) and (15), reach the same results. Therefore, the geometric parameters are involved in the equation, and d is the distance between the fracture and the upper surface of the specimen. It should be indicated that T 0 reflects the rock surface temperature so that it is the only parameter that can be directly obtained. Under actual conditions, it is a great challenge to determine the distance between the fracture and the surface with a fixed temperature. Accordingly, d can be defined as the vertical distance between the fracture and the rock surface at a fixed temperature. When d is infinitely small, T r ðxÞ%T 0 and Eq. (20) can be rewritten in the form below: In Eq. (21), parameters are simplified, thereby facilitating the description of local heat transfer to realize the heterogeneity of the seepage-heat transfer. This method can be used as a fully coupled seepage-heat transfer model II. In the next section, the two models will be verified and compared with the fixed heat transfer coefficient model.

Model validation
Zhao conducted a seepage-heat transfer experiment on granite with a single fissure (Zhao 1994). The entire experiment was performed in the Rock Thermodynamics Experimental System at the Imperial College of Technology. The test specimen was a cylindrical medium-coarse-grained granite with a length of 102 mm and a diameter of 51 mm. Figure 3 indicates that prior to the test, a single fissure was made by the tensile splitting. T 0 is the temperature of the outer surface of the rock sample, which is constant. Moreover, T in and T out denote the temperature of the inlet and outlet flow of the fracture, respectively. Seventy-eight various working conditions were considered in the experiment. Meanwhile, different initial temperatures varying from 90 to 140 C were considered for the rock. After reaching the steady-state condition, the water temperatures at the inlet and outlet ends were measured. Because of practical limitations, only the water temperature at the outlet was obtained and distribution of the water temperature along the surface of the fracture, and variations of the internal temperature of the rock were not be measured.
To evaluate the proposed model from viewpoint of the dynamic heat transfer coefficient, it is applied to simulate the experimental process. All simulations are carried out in COMSOL Multiphysics software. The seepage and the heat transfer modules of the software are combined and the seepage-heat transfer of the dynamic heat transfer coefficient is simulated. The model is meshed by triangular and quadrilateral. In this regard, 22465 triangular meshes, and 2860 quadrilateral meshes are used. The minimum and average unit mass of the model are 0.1357 and 0.821, respectively. The unit area ratio is 6.964 Â 10 À5 . Table 1 presents the physical parameters used in the simulation, which are mainly derived from published articles (Ma et al. 2020). The other geometric parameters are shown in Table 2.
In the established model, it is assumed that there is no phase change because of the temperature, variation. Therefore, 23 operating conditions with an initial temperature of 90 C are considered in the simulations. Then simulations are carried out using heat transfer coefficient models I and II. The obtained results in this regard are shown in Table 2. By comparison with the experimental data, the accuracy of the established model is verified. The obtained results from 23 test cases in Figure 4 show that there are differences in the fully coupled seepage heat transfer model, including two dynamic heat transfer coefficients. The outlet temperature calculated by model is overestimated, and the maximum relative deviation from the experimental data is 2.2%. This may be attributed to the parallel walls, where the water flowing through the fractures in the rock differs from that in the slab model. Moreover, the fractures in the rock may prevent the full development of speed and thermal layer. It is observed that the higher the fluid velocity, the smaller the deviation. The calculated results of the dynamic heat transfer coefficient using local heat transfer (model II) are in good agreement with the experiment and the relative deviations are less than 1%. Higher accuracy in this model mainly originates from smooth and parallel fractures in the simulation while the fractures in the experiment are morphological. Figure 5 shows the distribution of fluid velocity in a smooth fracture. There is a boundary layer in the fluid flow, and the velocity gradually increases from the boundary to the center. However, the actual velocity in the fracture changes in accordance with the crack morphology. Meanwhile, the influence of the crack position on the flow characteristics is not considered in model II. Figure 6 indicates that there are differences between the water temperature in the fractures obtained from dynamic and static heat exchanges. For comparison, the coupled model is calculated using the static heat transfer coefficient proposed by Zhao (Eq. (20)) and two dynamic heat transfer coefficients (model I and model II). In the simulations, the inlet temperature flow velocity and the fracture opening were set to 42 C, 10.63 mm/s and 19.17 lm, respectively. The water temperature in the     fracture calculated by the model I increases rapidly within 10 mm, and temperature does not change beyond 20 mm of the fracture. It is concluded that the heat transfer coefficient calculated by the boundary layer theory is higher than the experiment value so this theory overestimates variations in the water temperature. It should be indicated that this is in excellent consistency with the analytical results (Zhao and Geotechnics 2014). Figure 6 shows that the outlet temperatures obtained from the dynamic heat transfer coefficient used by model II and the fixed static heat transfer coefficient are the same, while the water temperature distributions along the fracture are quite different. It is found that the water temperature calculated by the model II increases rapidly within 10 mm, then increases relatively slowly and the temperature does not change beyond 50 mm of the fracture. Based on the obtained results from the fixed heat transfer coefficient, the water temperature increases steadily along the fracture. Considering the limitations of the experimental condition, there is no experimental on the distribution of the water temperature in fractures. Therefore, the water temperature distribution in the fractures calculated by different models should be regarded as the characteristics of the model used. There is a significant difference between the temperatures of the inlet water and the internal surface of the fracture, which produces a large temperature gradient. Moreover, the early period of the contact between low-temperature water and high-temperature rock fractures rapidly increases. Therefore, the water temperature changes obtained from model II are more in line with the cooling law. Compared to Eq. (20), the parameter calculations are considerably simplified in Eq. (21).

Parameter sensitivity analysis
In the previous section, it was verified that the dynamic heat transfer fully coupled seepage-heat transfer models are accurate, but the results of Eq. (21) are more consistent with the experimental results. In this section, it is intended to focus on the analysis of the established model. In this regard, Eq. (21) indicates that the heat transfer coefficient is related to the aperture and the injection velocity. To study the influence of parameters on the heat transfer, the initial temperature of the rock and the inlet water temperature are set to 90 and 42 C, respectively. Moreover, the time step length of the aperture and injection velocity are both 0.01. It should be indicated that the range of the aperture is from 0.1 to 3 mm, and the injection velocity varies from 0.01 to 0.1 m/s. Figure 7 shows the temperature distribution of the water flow along the length of the fracture under different apertures, and Figure 8 reflects the temperature distribution of the rock under different apertures. The following conclusions can be drawn from these two figures. The greater the aperture is, the smaller the increase in water temperature. When the aperture is 0.1 mm, the water temperature can increase to 87 C, but when the aperture increases to 3 mm, the water temperature can only increase to 58 C. As the aperture increases, more low-temperature water flows into the fracture to exchange heat with the rock, and the hydraulic gradient will remove more heat. As the aperture increases, the water temperature increases rapidly along the fractures and the range of change gradually decreases. The water temperature rises rapidly from the entrance of the fracture, but it rises slowly in the later stage.
With the injection of low-temperature fluid, convective heat transfer becomes the dominant heat transfer mechanism between the rock and the fracture, the fracture wall begins to cool down, and the cold front gradually expands into the rock. The rock temperature around the fracture entrance decreases rapidly, and the heat absorbed by the fluid gradually decreases. Therefore, a band-shaped low-temperature zone appears and extends toward the exit of the fracture. The shape of this zone depends on the variations of the aperture. As the aperture increases, the low-temperature area near the fracture increases and the amount of convective heat exchange between the fluid and the rock surface increases significantly. Accordingly, the rock temperature near the fracture is greatly reduced, and the maximum temperature drop is almost 30 C. As the aperture increases, the temperature gradient of the rock near the fracture and the water flow in the fracture increase slightly. Moreover, when the fracture is very narrow, variations of the temperature field by the seepage flow are low. When the aperture becomes greater, the influence of the thermal convection becomes significant. This is especially more pronounced in the middle area, where the temperature gradient gradually decreases. The temperature of the water outlet has a decreasing trend. However, as the aperture increases, the degree of influence on the temperature of the water outlet decreases. Figure 9 shows the temperature distribution of the water temperature along the length of the fracture under different injection velocities, and Figure 10 shows the distribution of the temperature field in the rock with different injection velocities. Figure 9 indicates that the faster the injection flow rate, the smaller the increase in the water temperature. For an injection velocity of 0.01 m/s, the water temperature can rise to 87 C, while when the injection velocity is set to 0.1 m/s, the corresponding water temperature can only rise to 56 C. Considering the velocity of the injected fluid, the contact time between the low-temperature water and high-temperature rock mass reduces. As the injection rate increases, the rising range of the water temperature in the early stage gradually decreases. In the case of a low flow rate, the water temperature at the outlet does not change significantly. This mainly originates from the small flow velocity, so that the convection heat transfer between water and rock is low. Accordingly, the heat exchange mechanism is not much different from the heat conduction inside the rock. Consequently, the influence of the real length of different fractures can be ignored. However, as the flow rate increases, the temperature of the fluid in different fractured channels at the outlet begins to vary. This may be attributed to the large fluid velocity so the convective heat transfer intensity is much greater than that of the conduction heat transfer. At this time, the true length of the fracture cannot be ignored. It should be indicated different lengths of the fluid circulation path result in different fluid temperature distributions. Figure 10 shows that as the injection velocity increases, the rock cooling area expands along the fracture direction. Moreover, as the water velocity increases, the corresponding temperature gradient of the rock near the fracture and that of water flow in the fracture decreases. Meanwhile, the degree of convective heat transfer between rocks and water on the fracture surface also increases significantly. Newtonian convective heat transfer equation indicates that when the temperature difference between the rock and fluid is constant, the amount of convective heat transfer is directly proportional to the size of the convective heat transfer coefficient h. On the other hand, Eq. (21) reveals that h is directly proportional to the fluid velocity. Therefore, as the water velocity increases, the corresponding heat convection between water and rock increases.
4.3. The influence of fracture seepage on the temperature field of tunnel surrounding rock 4.3.1. Calculation model In the calculation model, the surrounding rock of the tunnel is composed of a single fissure. In order to accurately study the influence of water seepage in the fracture on the temperature field of the surrounding rock of the tunnel, a single fissure rock mass in deep underground space is the research object. The calculation scope area is the tunnel surrounding rock with a single fissure, and the size is y Â z ¼ 40 m Â 25 m. The geometric model is shown in Figure 11.

Calculation parameters
The thermophysical parameters of the surrounding rock matrix of the simulated single-fracture tunnel are shown in Tables 3 and 4. The relevant parameters are from the literature (Kang et al. 2020).

Initial and boundary conditions
The initial temperature of the surrounding rock is 90 C. Accordingly, a constant temperature boundary with T ¼ 90 C is considered for the bottom boundary of the model. Moreover, the adiabatic condition is considered for other boundaries. Considering the influence of airflow, a constant temperature boundary with T ¼ 20 C is considered for the tunnel.
The initial condition of the seepage velocity in the fracture is set to 5 mm/s, and the non-slip condition is considered for the contact mode between the fracture and the surrounding rock. Then the velocity boundary at the entrance of the fracture is set. The coupled boundary condition of water-rock heat transfer in the fracture is the heat transfer coefficient at the junction of the fracture water and the surrounding rock matrix, which vary with temperature.

Mesh generation
The model is meshed by triangular and quadrilateral meshes, and the meshes are refined at the fracture. In this regard, 208128 triangular meshes and 22858 quadrilateral meshes are used in the model. The minimum and average unit mass of the model are 0.1355 and 0.8374, respectively. The unit area ratio is 2.156 Â 10 À6 .

Results and analysis
(1) Effect of the fracture aperture Figure 12 shows the temperature field distribution of the tunnel surrounding rock for the seepage velocity of 9 mm/s and different fracture widths. Considering the water flow in the fissures, the temperature field of the surrounding rock has different regional characteristics. For instance, the temperature of the surrounding rock of the tunnel is distributed along the fracture. It is observed that the temperature in the crack area is low and the temperature far from the crack area is high. This is because the flow of fissure water brings the heat energy of the surrounding rock downstream, causing the isotherms of the surrounding rock in Figure 12 to move to the right. When cold water is injected into the high-temperature rock, heat transfer occurs towards the fracture surface indicating that the seepage in the fracture directly controls the temperature distribution of the tunnel surrounding rock.
The temperature field between the fissure and the surrounding rock merges with the increase of the fracture width. As the width of the fracture increases, the temperature of the surrounding rock decreases and the fusion area, which is defined as the immediate cooling zone, expands. Then the cooling area and the surrounding high-temperature area form an asymmetric distribution of the overall surrounding rock temperature field. This Figure 11. Schematic diagram of computational model geometry. is because as the width of the fracture increases, the water content of the fracture in the surrounding rock increases, which causes the flow of the fracture water to take away more heat, thereby decreasing the temperature of the surrounding rock.
(2) Effect of the seepage velocity Figure 13 shows the temperature distribution of the tunnel surrounding rock or the crack width of 5 mm and different seepage velocities. It is observed that the temperature field of the surrounding rock of the tunnel forms different regions, resulting in an asymmetrical distribution of the temperature field of the surrounding rock. Figure  13 shows as the water velocity in the fracture increases, the asymmetric distribution range of the surrounding rock temperature field expands and the cooling zone moves to the right. It is concluded the seepage velocity of the fracture water has a significant effect on the distribution of the surrounding rock temperature. Moreover, it is observed that the low-temperature region of the fracture temperature field mainly moves to the right along the fracture. For the seepage velocity of 3 mm/s, an independent temperature field forms between the crack and the tunnel. When the seepage velocity increases to 6 mm/s, the fracture seepage temperature and the surrounding   Qu et al. (2017) 14 Zhao et al. (1992) 15  Lienhard and Catton (1986) 18 Lienhard and Catton (1986) 20 Zhao and Geotechnics (2014) rock temperature of the tunnel merge and form a joint zone. It is observed that when the seepage velocity reaches 9 mm/s, the cooling area is further enlarged and shifts to the seepage direction. This mainly originates from the increment of the water seepage velocity in the fracture, which requires more heat to dissipate, thereby decreasing the temperature of the surrounding rock. The faster the seepage velocity of the fracture water, the larger the cooling area of the temperature field, and the greater the disturbance to the temperature field of the surrounding rock.

Discussion on the efficiency of the convective heat transfer in fractures
The convective heat transfer efficiency of fractures is a core issue in the seepage-heat transfer. In engineering applications, the convective heat transfer efficiency of fractures can be improved by changing the hydraulic boundary conditions or improving the distribution of fractures. In this section, it is intended to calculate the convective heat transfer efficiency of the fractures under the dynamic heat transfer coefficient, and analyze the key factors affecting the convective heat transfer efficiency.

Variation of the overall convective heat transfer efficiency
Generally, the convective heat transfer efficiency of rocks can be described by the average convective heat transfer coefficient. In this regard, the average convective heat transfer coefficient of a fracture with a length of L and width of W can be expressed in the form below (Zhao 1994): where T r and T f are the average temperature of the rock fracture surface and the fluid, respectively. Moreover, A ¼ L Â W and q w denotes the heat flow in the fracture, which can be calculated by the fluid temperature difference and the flow rate between the inlet and outlet.
where q is the flow volume and q ¼ ubW: Figure 13. Influence of different seepage velocity in fracture on temperature field of tunnel surrounding rock.
Considering the dimensionless description of physical quantities, the convective heat transfer efficiency of the fracture can be described by the Nusselt number Nu.
Substituting Eqs. (22) and (23) into Eq. (24), the Nusselt number that characterizes the overall convective heat transfer efficiency of the fracture can be rewritten as follows: The higher the Nu number, the higher the convective heat transfer efficiency in the fracture. Figure 14 shows the distribution of Nu number for different injection velocities and fracture aperture conditions. It is observed that for a constant aperture, as the injection velocity increases, the corresponding Nu number of the fracture rises. The faster the injection velocity, the higher the variation rate of the Nu number. Meanwhile, when the injection velocity is constant, the Nu number of the fracture increases with the increase of the aperture. Combining the temperature distribution of the fracture along the flow direction in Figures 8 and 10 reveal that when a fluid with a lower temperature is injected into the fracture, the fluid temperature near the injection port changes significantly, while the outlet temperature almost remains constant. For a certain injection velocity, the larger the aperture, the larger the heat transfer flow.

Calculation and analysis of the local convective heat transfer efficiency of fractures
In the previous sections, it was found that the fluid temperature distribution has certain nonuniform characteristics. This feature can be used to study the distribution of the heat transfer efficiency of the fracture for uneven temperature distributions. In this regard, the fracture surface can be divided into several differential sections along the flow direction, corresponding to the average convective heat transfer coefficient. In this case, the local Nusselt number along the fracture can be defined as the following: where T 0 f ðxÞ is the derivative of the fluid temperature along the flow direction, indicating the local change in heat flow. Figure 15 shows the distribution of the convective heat transfer efficiency of the fracture along the flow direction for the aperture of 0.1 mm and different injection velocities varying from 0.01 m/s to 0.1 m/s. It is observed that when the dynamic heat transfer coefficient is used, the calculated Nu number along the fracture direction has a strong variability. Moreover, it is found that the Nu number, increases in a certain range close to the inlet, indicating that the convective heat transfer in this area is very significant. This phenomenon mainly originates from a large temperature gradient between the fluid and the rock when the fracture begins. As the heat exchange between rock and water progresses, the temperature difference between rock and water decreases so the corresponding Nu number decreases. Beyond the fracture direction of 60 mm and for the dynamic heat transfer coefficient, the Nu number drops below the static heat transfer coefficient Nu. Furthermore, the local Nu number of the fracture along the flow direction changes with the same law. Accordingly, it is concluded that the injection velocity has no significant effect on the local convective heat transfer efficiency of the fracture. This is consistent with previous investigations (Bai et al. 2017).

Conclusion
In the present study, a fully coupled model of seepage-heat transfer between a rock and flowing fluid in a fracture system was proposed based on the local dynamic heat transfer. Then a numerical simulation was carried out to evaluate the performance of the proposed model. The obtained results were compared with the experiment and it was found that the proposed model has excellent agreement with the experiment.
This model realizes the two-way coupling of seepage-heat transfer and, internally adjusts the heat transfer coefficient according to the dynamic changes of the system and covers the entire fracture. Controlling the aperture and water flow velocity are the most important factors affecting the temperature distribution of the fractured rock. It was found that with a wider aperture and faster flow velocity, the temperature of the entire rock decreases sharply.
The proposed model can reduce the number of parameters by locally calculating the heat transfer coefficient, thereby simplifying the calculations so that this method can be used in large-scale engineering applications. All parameters required in the model can be easily obtained through field observation, laboratory testing and numerical simulations.