Effects of fracture parameters and roughness on heat‐flow coupling in rock masses with two‐dimensional fracture networks

Discrete fracture network (DFN) is an effective means of describing the coupling of heat flow in an underground fractured rock mass. In this paper, an improved DFN is proposed by introducing the correlation function of fracture trace length and aperture, which is more consistent with the real fracture data. Next, based on the improved model, the influence of fracture roughness on the fluid flow and heat transmission was evaluated, and the relationship between the fracture aperture and the joint roughness coefficient (JRC) is established. Finally, based on the exponential function between confining pressure and aperture, the influence of confining pressure on the heat‐flow coupling process is considered in the improved model. Besides, the reliability of the model was verified by comparing with the analytical solution of the two‐dimensional single‐fracture heat‐flow coupling problem. The results show that under the same conditions, considering the correlation between the geometric parameters of the fracture, the seepage and heat transfer rates of the model increased, the outlet boundary flow reached the maximum value, and the average outlet temperature dropped rapidly. However, the fracture roughness reduces the outlet flow rate and the decline rate of average temperature. The confining pressure will lead to a decrease of about 3.5% in the outlet flow of the model, which is consistent with the seepage law in practical engineering. The model presented in this paper is an effective supplement to the two‐dimensional fracture network heat‐flow coupling model and can provide a theoretical basis and a numerical calculation tool for related underground rock mass engineering.


| 1217
HUANG et Al fractures. Fractures are considered as the main route for fluid flow and solute transport. 5,6 In addition, in the process of geothermal exploitation, a large number of fractures are generated in the underground rock mass through hydraulic fracturing. The flow of fluid in the fractures changes the temperature of the rock mass. 7 Therefore, the study of heat transfer and percolation mechanism in rock mass fractures is an important topic that scholars in the field of hot dry rock analysis and rock mass hydraulics have focused on in the recent years.
To effectively simulate the heat transfer and flow of underground fractured rock masses, many scholars have developed and established a variety of numerical simulation tools and theoretical models. [8][9][10] Chen et al 11 established a porous media finite element model for multiphase flow heat transfer and stress-strain coupling according to the law of momentum mass and energy conservation in continuum mechanics. Tong et al 12 presented the development of an effective thermal conductivity model for the simulation of thermohydro-mechanical processes in geological porous media. The derivation of a new effective thermal conductivity model in closed form is then presented. Baoli et al 13 proposed the finite-size particle method (FSP) and double-distributionfunction (DDF) coupling model of the lattice Boltzmann method for fluid flow and heat transfer in porous media. The above research can simulate the heat-flow coupling process of fractured rock mass at a macroscopic level. However, the information about the fracture structure was ignored, and the distribution of fractures in the rock mass could not be studied. Meanwhile, the influence of the variation of the geometric parameters such as fracture aperture, trace length, and roughness on the seepage field and thermal field is not considered, so the calculation results are often idealized. Therefore, it is necessary to characterize the fracture and analyze the effect of the fracture geometry parameters on the heat-flow coupling.
In the 1970s, the discrete fracture network (DFN) model was introduced to characterize the fractures in rock mass, and it has been widely used because of its intuitive effectiveness. 14, 15 Chen et al 16 18 focused on the construction of a rough-walled discrete fracture model to simulate mass and heat transfer in a geothermal reservoir by using water and CO 2 as the working fluid. However, although the DFN model is widely used to simulate the thermal-hydraulic-mechanicalchemical coupling process of rock mass containing a large number of fractures, in most previous DFN models, the geometric parameters of fractures, such as the fracture location, trace length, direction, and hydraulic aperture, are generally considered to be irrelevant.
The above studies also did not consider the correlation between fracture trace length and aperture when simulating random fracture networks. Some researchers inferred that there was a certain correlation between the trace length and aperture of the fracture through the observation of field mapping results. [19][20][21] This correlation has a great influence on the formation of the fracture network and heat-flow analysis, and even leads to a large error in the analysis results. In addition, the fracture in the numerical model is usually simplified as an idealized smooth fracture without considering the effects of the fracture roughness. 17 As shown in Figure 1, the roughness has a significant influence on seepage. To date, the effect of fracture roughness on the heat-flow coupling process needs further study. 22,23 In view of this, an improved fracture network model is established by introducing the correlation function between fracture trace length and aperture; the influence of fracture roughness and confining pressure on the fluid flow and heat transmission are evaluated. The simulation of a more practical fracture field can supplement the development of a stochastic fracture network model, which has a F I G U R E 1 A, Rough discrete fracture network (left, varying apertures in each single fracture); (B) idealized discrete fracture network (right, smooth, and parallel-walled fractures). Line thickness represents the magnitudes of fracture apertures guiding significance for underground nuclear waste repository, carbon dioxide sequestration, and enhanced geothermal systems.

| Assumptions of fracture geometry parameters
The locations of fractures in DFN simulations are usually assumed to follow Poisson's distribution. 24 In order to reflect the location of fractures, the center point of fractures obeying Poisson distribution is used to represent. The Poisson process is realized by random numbers generated by recursive algorithm. The recursive formula is as follows: where R i is the uniform number of random distribution between [0,1]; INT(x) is the integral function. This paper takes the x integer part of the internal function. The initial value R 0 is generated by the linear congruence method. Figure 2 shows the location coordinates of the fracture center in the DFN numerical model with a study area of 30 × 30 m. The fracture density in the figure is 0.45/m 2 .
In the field outcrop survey, the literature 25 gives the cumulative probability density function of the power-law distribution of fracture trace length by analyzing the geometric characteristics of some relatively large fractures: where L is the fracture trace length; L max and L min are the maximum and minimum fracture trace lengths; the value range of L is 0.5-250 m; F is a random uniform distribution with interval of [0,1] 25 ; D is the fitting fractal dimension 26,27 ; and its value range is 2.2 ± 0.2. Figure 3 shows the cumulative probability curve of fracture trace length with different fractal dimensions D.
According to literature, 24,27 fracture orientation follows Fisher's distribution, which is widely used in numerical modeling of discrete fracture networks. The cumulative probability density function of deviation angle of average orientation is given by the following formula: where K is Fisher's constant, and F is a random uniform distribution between [0, 1]. It can be seen from Figure 4 that the deviation angle ϴ decreases with the increase of K. When K = 10, nearly 75% of the fractures of the deviation angle ϴ are less than 0.5 rad, indicating that the overall variation of the fractures in the region is small. However, when K = 1 and K = 5, the deviation angle ϴ is generally large and some even reach 3 rad, which is contrary to the original orientation of the fractures.
The fracture aperture can be indirectly calculated from the laboratory shear-flow test by assuming that the fluid movement of a single fracture follows the cubic law, and its calculation formula is as follows 27 : where h is the fracture aperture, w (m) is the flow passage width, V (m 2 /s) is the fluid motion viscosity coefficient, Q (m 3 /s) is the flow rate, i is the hydraulic gradient, and g (m 2 /s) is the gravitational acceleration. To simplify the calculation, the fracture aperture is taken as a constant value. According to the standard hydraulic test results, 28,29 the fracture aperture obeys the lognormal distribution, and its density function is as follows: where a and b are the first and second moments of lognormal distribution. In the numerical simulation of DFN using the fracture geometry data obtained from field outcrop survey, the distribution function of truncated fracture trace length and aperture is unavoidable. In this study, the truncation distribution function of fracture aperture is as follows 30,31 : where h a and h b are the upper and lower limits of the fracture aperture, respectively. The range of fracture aperture in this paper is 1-200 μm. When the fracture aperture is less than h a and greater than h b , the value of the truncated distribution function is 0. The truncated cumulative distribution function of lognormal distribution is as follows: After algebraic operation, the above equation can be written as: is an error function. From Equation (4) to Equation (8), the fracture aperture distribution function can be written: where erfinv(x) is the inverse error function, and F is a random number between [0,1].

| Relationship between fracture trace length and aperture
For simplicity, fracture trace length and aperture are usually assumed to be unrelated. However, when some scholars studied the fracture modes in the process of clay tensile deformation, they found correlation characteristics between fracture trace length and aperture, [19][20][21] and proposed a power-law correlation function between fracture trace length and aperture: where β is the characteristic index, and its value is between 0.5 and 2. Upon studying the correlation between fracture trace length and aperture, 30 it was pointed out that with the increase of fracture trace length and displacement, the average fracture aperture tends to increase, and the truncated cumulative probability density function of fracture trace length can be defined as 31 : where L max and L min are the maximum and minimum fracture trace lengths, respectively. By substituting Equation (11) into Equation (9), a correlation function based on the logarithmic normal distribution function of fracture aperture and the power-law distribution function of fracture trace length can be obtained: where D is the fractal dimension of the fracture, h a and h b are the upper and lower limits of the fracture aperture size, respectively, and h is the fracture aperture. The initial value of h is given by Equation (4). Therefore, a and b can be directly used as input parameters, and then, the fracture aperture can then be calculated by Equation (12). A large amount of experimental data shows that the fluid flow in rock mass fractures is greatly affected by stress, which is mainly because the geometrical parameters of fractures change with the stress field. 32, 33 Sun et al 34 proposed that the change in hydraulic aperture e h with confining stress σ n can be represented by the following exponential function: where e 0 is the initial hydraulic aperture, and α is the fracture deformation parameter. By introducing the above formula into the numerical simulation of the thermal flow coupling of fractured rock mass, the fracture aperture is made to change with time under the confining stress, in order to study the influence of confining pressure on the heat-flow coupling process of fractured rock mass.

| Fracture roughness
Since Barton's proposed joint roughness coefficient (JRC) in the empirical criterion of shear strength of joints, 35 most studies have been devoted to accurately and objectively measuring JRC values of fracture surfaces. However, little attention has been paid to the JRC distribution of fractures in the study area, which is mainly due to the lack of field data. Luo et al 17 observed that JRC was fitted to some common probabilistic distributions, that is, normal, lognormal, Weibull, and Gamma distributions, through the study of two databases. Based on the above studies, this study only considers the effect of fracture surface roughness on fluid flow and heat transfer in the case of a lognormal distribution.

| Determination of fracture hydraulic aperture
In this study, the influence of different factors on the heatflow coupling process of a fracture network is reflected by the fracture aperture. Therefore, it is necessary to determine the fracture hydraulic aperture applied in the heat-flow coupling process before modeling. In this study, Equation (12) is the mechanical aperture obtained according to the mechanical tests, and the mechanical aperture of rock mass fractures is usually larger than the hydraulic aperture, which is mainly due to the effect of fracture surface roughness. 36 Therefore, it is necessary to introduce an empirical formula between the hydraulic aperture e h and the mechanical aperture e m to represent the influence of the fracture surface roughness. Zhao and Li 37 review the main empirical models between hydraulic aperture e h and mechanical aperture e m , and consider the parameters σ/E, JRC, and Z 2 , as shown in Table 1. The empirical model between the hydraulic aperture e h and the mechanical aperture e m proposed by Li and Jiang 38 is considered in this paper, and its empirical formula is as follows: Grasselli et al 39 studied the shear strength of rock joints, proposed the concept of JRC by comparing the shape of natural rock joints, and defined the standard section curves of joint roughness, as shown in Figure 5. The empirical expression is as follows: Re > 1

Researchers e/E
Patir and Cheng Zimmerman and Bodvarsson where Z 2 is the root mean square of the slope of the fracture plane profile. Through Equations (14) and (15), the relationship between the roughness coefficient JRC and the ratio of hydraulic aperture e h and mechanical aperture e m was established. The mechanical aperture e m can be calculated using Equation (12). Thus, the hydraulic aperture can be obtained by considering the relationship between fracture trace length and aperture, confining pressure and the influence of fracture roughness and other factors. The distribution function of roughness coefficient JRC is introduced to evaluate the influence of roughness on rock thermalflow coupling in this paper, which is more representative.

| Governing equations for fluid flow and heat transfer
A fractured rock mass is a binary medium composed of the rock mass matrix and the discrete fractured rock mass. In this paper, COMSOL Multiphysics was used to establish a DFN model to simulate the fluid flow and heat transfer in rock mass fractures.

| Governing equation of rock matrix
The mass conservation equation of the fluid flow in the rock mass matrix is 40 : By substituting the above equation into the continuity equation, the seepage field equation with pressure as the dependent variable can be obtained as 40 : Assuming that the water flow in the rock mass matrix obeys Darcy's law, the water flow rate can be expressed as: where t (s) represents time, ρ (kg/m 3 ) is the density of water, S m (1/Pa) is the water storage coefficient of the rock mass matrix, u (m/s) is the flow rate of the fluid in the rock mass matrix, ε p is the porosity of the rock matrix, κ (m 2 ) is the permeability of rock matrix, μ (Pa*s) is the viscosity of water, and p (Pa) is the pressure.
The heat field in the rock mass matrix can be characterized using the energy conservation Equation 40 : where T (K) represents the temperature; k (W/(m* K)) is the thermal conductivity of the fluid; k eff (W/(m*K)) is the thermal conductivity of the rock mass matrix; c p (J/(kg*K)) is the specific heat capacity of the fluid; k p (W/(m*K)) is the thermal conductivity of the rock mass matrix; ρ p (kg/m 3 ) is the matrix density of the rock mass; c p.p (J/(kg*K)) is the specific heat capacity of the rock mass matrix; and Q md (W/m 3 )is the heat source.

| Fracture control equation
The governing equation of the seepage field in fractures is similar to that of the rock mass matrix, and the continuity equation of fractures can be described as follows 40 : Assuming that the water flow in the fracture also obeys Darcy's law, the water flow velocity is as follows: are the pressures on both sides of the fracture along the normal direction, respectively; and the fracture has a small aperture. Therefore, it is assumed that the pressure on both sides of the normal fracture is similar, that is: where p f is the pressure inside the fracture. Fracture heat transfer balance equation is given by 16 : In Equation (26), the first term is the change in energy in the fracture, the second term is the increase of the fracture energy by convection, and the third term is the increase in energy in the thermal conductivity fracture. In this equation, f e up and f e bottom represent the energy transferred from the fracture method to the fracture on each side, respectively.

| Establishment of the model
Geothermal reservoirs are buried deep underground, and the detailed geometry of fractures is difficult to observe directly, which poses a challenge to numerical modeling. Therefore, the fracture network generated randomly in this study was used to characterize the geothermal reservoirs and simulate the heatflow interaction process in the geothermal heat production process, with a model size of 100 × 100 m. The parameters of the discrete fracture network generated using the method introduced in Section 2 are shown in Figure 6. It can be seen that the fracture strike distribution in this region is random and the fracture length is mainly between 2.5 and 3.5 m. Figure 6C shows the flow path after eliminating isolated or closed fractures, and the fracture endpoint x 0 is the fracture outlet of the flowing water. At the same time, to maintain the simplicity and repeatability of the model, certain initial conditions and boundary conditions are assumed: 1. Seepage field: Preset pressure boundary conditions were adopted to ensure water circulation in the system. It is assumed that water flow is injected from the left boundary of the model at x = 0 m with a constant water head of 1000 m and the right boundary at x = 100 m with a constant water pressure of 0 MPa. 2. Thermal field: Let us assume that the initial temperature of the fractured rock mass be 85°C. The x = 0 m on the left-hand side is the waterflood boundary, and the initial temperature at the inlet is T in = 20°C. On the right-hand side, x = 100 m is the outlet boundaries, and the other boundaries are the thermal insulation boundary. The total calculation time was 30 d, and the step length is 1 d using the transient simulation solution.
By discretizing the rock mass matrix and fracture into finite triangle elements and line elements in COMSOL, respectively, and combining Sections 3.1, and 3.2, the following four two-dimensional rock mass models with different fracture apertures are established: A The apertures of all the fractures are 0.065 m, using this model as the reference group.  (14) and (15) are introduced to consider the effect of fracture surface roughness in the process of heat-flow coupling. The average fracture aperture was 0.0754 m. D Based on the model of C), Equation (13) is introduced for the change in hydraulic aperture e h with confining stress σ n . It is assumed that the hydraulic aperture also changes with time. To facilitate the modeling, MATLAB software was used in this paper to calculate the numerical value of the fracture aperture varying with confining pressure within 30 d. The value was assigned to the fracture in the subsequent heat-flow coupling calculation process, so that the influence of confining pressure on the heat-flow coupling process can be evaluated. The parameters of the numerical model are presented in Table 2. Figure 7 shows the curve of fracture aperture varying with time during the model D calculation.

| Verification of the model
The model needs to be verified to ensure the reliability of the calculated results before the simulation. In this study, the numerical model was verified by comparing it with the analytical solution of the two-dimensional single-fracture heat-flow coupling problem. This method has been widely used. 40,41 As shown in Figure 8, a horizontal fracture with a constant aperture divides a homogeneous, isotropic, and impermeable rock into two parts. At the same time, a Cartesian coordinate system was set to make the origin coincide with the entrance.
In the single-fracture rock mass model, the left end of the fracture is the inlet and the right end is the outlet. Assuming that the width of the rock extends to ±∞ in the y direction, the length of the rock matrix and the fracture length extend to infinity in the X direction, forming a semi-infinite coordinate system. Assuming that the initial temperature of the singlefractured rock mass model is T 0 , in the process of heat-flow simulation, water is injected into the model from the left end at a constant velocity u f and a constant temperature T in , and then flows out from the right end of the fracture. The fracture aperture is much smaller than the fracture length, so the temperature change of the fracture along the y-axis was not considered. Barends 42 provided analytical solutions for the temperature changes of fractures at different locations in a two-dimensional single-fracture model: where T f is the internal temperature of the fracture; ρ is the density; c is the specific heat capacity; λ is the thermal conductivity; the subscripts s and f are the rock mass matrix and water, respectively; erfc is the complementary function of error; and U is the unit step function.
To make a good comparison with the above analytical solutions, a 100 × 100 m rectangular study area was used to characterize the rock mass matrix of a single-fracture model in the calculation of the numerical model. Figure 9 shows the numerical simulation of the fracture network model. The initial temperature of the rock mass matrix was 85°C, the fracture inlet temperature was 20°C, and material parameters are listed in Table 3. Figure 10 shows the comparison between the analytical and numerical solutions of the single-fracture model. Figure 10A shows the temperature variation with time for three different locations in the fracture, x = 40, 60, and 80 m, and (28)

Time (days)
the temperature is the centigrade scale. Figure 10B shows the temperature variation with position at three different times in the fracture, at time t = 10, 50, and 100 d, where d is used to represent the time in days. It can be seen from the figure that the numerical solution is basically the same as the analytical solution except for a slight deviation. This bias may be due to the fact that the analytical solutions are based on the assumption of semi-infinite systems, whereas the numerical solutions are derived from finite geometric models. Therefore, the numerical model calculation of the heat-flow coupling process is reliable.

| Results of the hydraulic process
The percolation process of underground fractured rock mass is very important for geothermal exploitation and nuclear waste burial, and the flow velocity of the fluid in an underground fractured rock mass may directly affect the heat transfer process and nuclear waste diffusion process. Figure 11 shows the calculation results of the seepage process of the four models. Figure 11A shows the change for outlet boundary flow rate Q with time in different models. It can be seen that the flow rate at the outlet boundary of model A is the smallest, which is caused by the small aperture of the fracture. The flow rate at the outlet boundary of model B is the largest, mainly because the assumed uniform fracture aperture is less than the fracture aperture calculated by Equation (12), and its average aperture is much larger than that of the other models. Therefore, the results of the seepage calculation without considering the correlation between trace length and aperture often have a large error. Based on model B, model C considers the influence of the normal distribution roughness coefficient JRC on the seepage flow. It can be seen that after considering the fracture surface roughness, the outlet boundary flow rate relatively decreases. This also indicates that the mechanical aperture is often larger than the hydraulic aperture. This is because the surface roughness of the fracture leads to a tortuous effect and eddy current in the flow process. Model D assumes that the fracture aperture changes with time based on model C. It can be seen that the flow rate at the outlet boundary decreases with time, and finally stabilizes at a smaller value. This phenomenon is often consistent with actual seepage in a fractured rock mass. Figure 11B shows the change in the flow rate with time at fracture endpoint x 0 . It can be seen that with the decrease in the fracture aperture, the flow rate in the fracture also decreases. This phenomenon is consistent with the flow law and the seepage as shown in Figure 11A. Figure 11C shows the change in pressure with position on the path y = 50 m at t = 30 d. The hydraulic pressure of model A on the path is the largest, because the fracture aperture of model A is the smallest. According to Darcy's law, the water pressure drops slowly along the flow direction in model A. At the same time, it can be seen that there is little difference in pressure change among models B, C, and D. This shows that the influence of fracture surface roughness on pressure or the influence of the fracture aperture variable on pressure has little influence on numerical model calculation results in this study. This is mainly due to two reasons. First, the seepage pressure was mainly affected by the fracture aperture and permeability. Second, in this study, the fracture aperture showed little change, all within the same order of magnitude, and had limited influence on the seepage process. Similarly, it can be seen in Figure 11D that the pressure changes with position on the path x = 50 m when t = 30 d.  Figure 12A shows the seepage pressure simulation results of model D at t = 30 d. It can be seen that the seepage pressure in the study area is affected by the fracture location. Figure 12B shows the velocity distribution in the study area, and it can be seen that the flow velocity increases significantly between the two fracture endpoints. This is because the permeability of the fracture is much greater than that of the rock matrix. Figure 13 shows the  Figure 13D shows the change of pressure with time at three different positions (x = 25 m, x = 50 m, and x = 75 m) in the discrete fracture network. It can be seen from Figure 13A that the head pressure at t = 30 d is lower than that at t = 1 d. At the same time, it can be seen from Figure 13B that the head pressure changes are different due to the inconsistent distribution of fractures in different paths. It is found that the discrete structure of the fracture network can control the seepage and thermal fields significantly. Figure 14 shows the temperature distributions of different fracture network models at different times. By transverse comparison, it can be seen that at the beginning of the simulation, owing to the injection of low-temperature water, the rock matrix and fracture underwent heat exchange with the flowing fluid, and the fluid temperature immediately increased. The temperature of the rock matrix and fracture also decreased gradually, forming a low-temperature zone. As more heat is extracted from the reservoir, the cryogenic zone expands until a global cryogenic zone is formed. At the same time because the permeability of the fracture is much larger than that of the rock mass matrix, the decrease in the temperature of the rock mass matrix is slower than that of the fracture. It can be observed that the temperature drops more rapidly near some fracture concentration areas. These connected fractures form the main flow channel of the circulating fluid. In connected fractures, the effect of thermal convection is very strong owing to the high flow velocity. The interaction of the heat transfer process is consistent with the results in the literature, 40,41 and it shows that the fracture distribution has a significant influence on the heat transfer of the DFN model. Therefore, the influence of fracture distribution and geometric parameters must be considered in both heat transfer and seepage analyses.

| Results of the thermal process
Comparing the temperature distributions of the four different models simultaneously, it can be seen that the temperature distribution of model A with the uniform fracture aperture differs little from that of model D. The fracture aperture of model A is smaller than the average fracture aperture of model D. This shows that roughness and confining pressure have a significant influence on the heat transfer process. Model B assigns the aperture value of each fracture separately, and it can be seen that it has the highest heat transfer rate. It can be seen from the temperature distribution diagram of model C that the effect of fracture surface roughness on heat transfer cannot be ignored. The simulation results of model D show that the heat transfer rate decreases with a decrease in the fracture aperture. Figure 15 shows the temperature variation with position at three different times in the DFN. It can be seen that the temperature distribution is consistent with that in Figure 14. It is worth noting that when t = 30 d, the temperatures of models A and D are higher than that of the other models, mainly because of the small fracture aperture of model A and the influence of the change of fracture aperture of model D. Figure 16 shows the changes of the average outlet temperature and the temperature at the fracture end x 0 , respectively. The average outlet temperature of model A is slightly higher than that of the other models, mainly because of its low fracture aperture. The average outlet temperature of model C is slightly higher than that of model B. It reflects that the roughness of fracture surface hinders the heat transfer process. This is also consistent with the law that the mechanical aperture is larger than the hydraulic aperture. 37 From the temperature distribution curve of model D, it can be seen that with the increase of time, the fracture aperture decreases, which leads to the decrease in the heat transfer rate of the fracture.

| CONCLUSION
Four different heat-flow coupling models were established by programming the grid generation program in MATLAB. The effects of fracture parameters, roughness, and confining pressure on the thermal flow coupling of fractured rock masses are discussed. The conclusions are as follows: The discrete structure of the fracture network has an obvious control effect on the seepage and thermal fields. It is necessary to further elucidate the influence of the fracture geometry parameters on the coupled thermal flow analysis of the rock mass. Under the same conditions, considering the correlation between the geometric parameters of the fracture, the seepage and heat transfer rates of the model increased, and the outlet boundary flow reached the maximum value. The heat transfer velocity is faster than the traditional model, which makes the heat-flow coupling process more conform to engineering practice. The improved fracture network model can complement the development of random fracture network model. The fracture surface roughness can reduce the outlet flow. At the same time in the heat transfer process, the average outlet temperature decreased relatively slowly. The confining pressure can reduce the outlet flow, which is consistent with the seepage law under creep conditions. This indicates that the confining pressure will slow down the heat transfer process of the model, which cannot be ignored in practical engineering.
This two-dimensional model provides a means for in-depth study of the heat flow of deep geothermal reservoirs and is helpful for improving the design of underground nuclear waste storage and geothermal reservoirs. In future works, it will be necessary to study a three-dimensional numerical model.