## Numerical Simulations  
### Mathematical and Physical Models  
A three-dimensional single fracture model was established to study the influence of flow direction on the heat transfer characteristics in a granite single fracture. The geometry profile of the fracture was constructed using a random method (Figure 1a), while the rough fracture wall was built by stretching in the vertical direction (Figure 1b,c). The surface flow of the fracture surface was assumed to have four different directions, with the 0\({}^{\circ}\) direction taken as the \(z\)-axis direction (Figure 1b). The angle between the surface flow direction and the positive direction of the \(z\)-axis is denoted as \(\alpha\). The surface flow direction had a gradient of 30\({}^{\circ}\) and was rotated counterclockwise around the \(z\)-axis by 30\({}^{\circ}\), 60\({}^{\circ}\), and 90\({}^{\circ}\) to obtain the fracture wall. The rough single-fracture numerical model shown in Figure 2 was obtained by embedding these fracture walls into a cylindrical model with a size of 100 \(\times\) 50 mm (length \(\times\) diameter). The vertical displacement of the fracture in the model was determined as 0.3 mm by translation method. Table 1 lists the geometric parameters of the four fracture models. The surface area and inlet area  
Figure 1: Fracture surface and surface flow direction distribution. (**a**) 2D random fracture profile; (**b**) pseudo-three-dimensional fractured rough surface stretched towards the page direction. (**c**) 3D view of pseudo-three-dimensional fracture rough surface. deviations of the four fractures were 1.87% and 0.33%, respectively. Therefore, at the scale of this study, the effects on the simulation results due to differences in surface area and inlet area are negligible. For natural fracture, the Reynolds number is much smaller than the critical Reynolds number (Re \(<\) 1800 [35]) due to its small aperture and low flow rate of underground fluid. Thus, the fluid flow in a granite fracture remains laminar. Convective heat transfer in three-dimensional rock fractures was modeled using FLUENT 6.3 [40] by employing a finite volume method (FVM) to solve the Navier-Stokes equations. The governing equations for the conservation of momentum, mass, and energy in the fracture are described as follows:  
Continuity equation:  
\[\frac{\partial\rho_{w}}{\partial t}+\
abla\cdot(\rho_{w}\mathbf{V})=0 \tag{1}\]  
Momentum equation:  
\[\frac{\partial(\rho_{w}\mathbf{V})}{\partial t}+\
abla(\rho_{w}\mathbf{V})=-\
abla P+ \
abla\cdot\left(\mu\left[\left(\
abla\mathbf{V}+\
abla\mathbf{V}^{T}\right)-\frac{2} {3}\
abla\cdot\mathbf{V}I\right]\right) \tag{2}\]  
Energy equation:  
\[\frac{\partial}{\partial t}(\rho_{w}E)+\
abla\cdot(\mathbf{V}(\rho_{w}E+p))=\
abla \cdot\left(\lambda_{w}\
abla T+\mu\left[\left(\
abla\mathbf{V}+\
abla\mathbf{V}^{T} \right)-\frac{2}{3}\
abla\cdot\mathbf{V}I\right]\cdot\mathbf{V}\right) \tag{3}\] \begin{table}
\begin{tabular}{c c c c c c} \hline \hline
**Model Case** & **0\({}^{\circ}\)** & **30\({}^{\circ}\)** & **60\({}^{\circ}\)** & **90\({}^{\circ}\)** & **Maximum Deviation** \\ \hline Surface area & & & & & \\ (mm\({}^{2}\)) & 5588.95 & 5620.33 & 5658.53 & 5553.84 & 1.87\% \\ Inlet area & & & & & \\ (mm\({}^{2}\)) & 14.9883 & 14.9391 & 14.9405 & 14.9425 & 0.33\% \\ Volume & & & & & \\ (mm\({}^{3}\)) & 1496.28 & 1495.96 & 1494.91 & 1494.25 & 0.14\% \\ Vertical offset & & & & & \\ (mm) & & & & & \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Geometric parameters of the fracture. Figure 2: Fracture seepage heat transfer model. where \(\rho_{w}\) (kg/m\({}^{3}\)) and \(\rho_{r}\) (kg/m\({}^{3}\)) are the densities of water and rock, respectively; \(V\) (m/s) is the velocity vector; \(p\) (Pa) is the pressure; \(I\) is the unit tensor; \(E\) is the total energy; \(\lambda_{w}\) (W/(m\(\cdot\)\({}^{\circ}\)C)) is the fluid thermal conductivity; \(\mu\) (Pa s) is the fluid viscosity. In EGS frameworks, water is typically used as a heat exchange medium.