<figure><img src="image_1.png" /></figure>

HAL
open science

# Heat transport by flow through rough rock fractures: a numerical investigation

Maria Klepikova, Yves Méheust, Clément Roques, Niklas Linde

► To cite this version:

Maria Klepikova, Yves Méheust, Clément Roques, Niklas Linde. Heat transport by flow through rough rock fractures: a numerical investigation. Advances in Water Resources, 2021, 156, pp.104042. 10.1016/j.advwatres.2021.104042. insu-03348788

HAL Id: insu-03348788

https://insu.hal.science/insu-03348788

Submitted on 20 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Journal Pre-proof

Heat transport by flow through rough rock fractures: a numerical
investigation

Maria Klepikova, Yves Méheust, Clément Roques, Niklas Linde

PII: S0309-1708(21)00196-2

DOI: https://doi.org/10.1016/j.advwatres.2021.104042

Reference: ADWR 104042

To appear in: *Advances in Water Resources*

Received date: 19 July 2021

Revised date: 14 September 2021

Accepted date: 14 September 2021

Please cite this article as: Maria Klepikova, Yves Méheust, Clément Roques, Niklas Linde, Heat transport by flow through rough rock fractures: a numerical investigation, *Advances in Water Resources* (2021), doi: https://doi.org/10.1016/j.advwatres.2021.104042

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2021 Published by Elsevier Ltd.

<figure><img src="image_3.png" /></figure>
1 highlights

2
3 We analyse advective and conductive heat transport in rock fractures, coupled to heat conduction in the matrix

4
5 Flow channeling strongly impacts heat transport through fractures and fracture-matrix heat exchange

6
7 The thermal behavior at the fracture scale can be predicted from the fracture's effective hydraulic transmissivity
# Heat transport by flow through rough rock fractures: a numerical investigation

Maria Klepikova<sup>a,b</sup>, Yves Méheust<sup>b</sup>, Clément Roques<sup>b</sup>, Niklas Linde<sup>a</sup>

<sup>a</sup>organization=University of Lausanne, Institute of Earth Sciences, city=Lausanne, postcode=1015, country=Switzerland

<sup>b</sup>organization=Univ. Rennes, CNRS, Géosciences Rennes, UMR 6118, city=Rennes, postcode=35000, country=France

## Abstract

Fracture surface topography exhibits long-range spatial correlations resulting in a heterogeneous aperture field. This leads to the formation, within fracture planes, of preferential flow channels controlling flow and transport processes. By means of a 3-D heat transport model coupled with a 2-D fracture flow model based on the lubrication approximation (i.e., local cubic law), we investigate how the statistical parameters determining spatial aperture variations in individual fractures control the heat exchange at the fluid/rock interface and heat transport by flow. Ensemble statistics over fracture realizations provide insights into the main hydraulic and geometrical parameters controlling the hydraulic and thermal behaviour of rough fractures. Similarly to the rough fracture's hydraulic behaviour, we find that its heat transport behaviour deviates from the conventional parallel plate fracture model with increasing fracture closure and/or decreasing correlation length. We demonstrate that the advancement of the thermal front is typically slower in rough fractures compared to smooth fractures having the same mechanical aperture. In contrast with previous studies that neglect temporal and spatial temperature variations in the rock matrix, we find that the thermal behavior of a rough-walled fracture can, under field-relevant conditions, be predicted from a parallel plate model with an aperture equal to the rough fracture's effective hydraulic aperture. This greatly simplifies the prediction of possible reservoir thermal behavior when using field measurable quantities and hydrological modeling.

Keywords: Fracture, Roughness, Heat exchange, Flow channeling
1. Introduction

Heat transport in fractured media is often considered in hydrogeological studies, for instance, when inferring hydraulic parameters by fitting heat transfer equations to thermal data. Heat carried by groundwater serves as a tracer that can be used to quantify flow through fractures [Ge, 1998; Read et al., 2013], to characterize fracture network connectivity [Silliman and Robinson, 1989; Klepikova et al., 2011, 2014] and to constrain regional scale flow patterns [Anderson, 2005; Saar, 2010]. Understanding heat transport in fractured media is a prerequisite for studying hydrothermal flows [Fairley, 2009; Malkovsky and Magri, 2016]. Moreover, characterizing heat transport processes in the subsurface is essential for numerous industrial applications. For instance, heat transfer is critical when assessing heat storage in the ground [Lanahan and Tabares-Velasco, 2017; de La Bernardie et al., 2019] and near-field thermal effects in the context of radioactive waste disposal [Zhang et al., 2017]. Knowledge of thermal transport is also necessary to maximize the efficiency and sustainability of geothermal systems [Kolditz and Clauser, 1998; Martinez et al., 2014; Shortall et al., 2015; Guo et al., 2016; Vik et al., 2018; Patterson and Driesner, 2020].

Heat transport in fractured media has predominantly been addressed using simplified conceptual fracture models. For example, most fracture network models used for geothermal investigations assume a 1-D linear flow geometry [Pruess and Doughty, 2010], or consider fractures as parallel-plate systems with a constant aperture [Gringarten et al., 1975; Kolditz and Clauser, 1998; Kocabas, 2005; Jung and Pruess, 2012; Zhou et al., 2017; Vik et al., 2018]. Hydrothermal studies generally represent faults as tabular bodies of internally homogeneous properties or as 2-D discontinuities that juxtapose hydrogeologic units of differing properties [e.g. Malkovsky and Magri, 2016]. While it is common practice to neglect fracture heterogeneity, the consequences of such simplifications in terms of predictability remain poorly understood [Klepikova et al., 2016; de La Bernardie et al., 2019].

The fracturing process itself, as well as post-fracturing processes such as geological stress and strain, chemical dissolution, precipitation and erosion, may result in complex fracture-wall surface geometries. At the scale of a single fracture, fracture wall roughness exhibits long range spatial correlations that induce a heterogeneous aperture field [Brown, 1987; Johns et al., 1993; Candela et al., 2012], thus, promoting strongly heterogeneous flow path distributions and preferential flow channels within fracture planes [e.g. Tsang
54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

and Tsang, 1987; Méheust and Schmittbuhl, 2000; Méheust and Schmittbuhl, 2001]. Numerous studies have shown that fracture roughness has a remarkably strong impact on fluid flow and, as a consequence, on solute and particle transport through single fractures. Studies of flow and solute transport in rough-walled fractures include both theoretical and numerical studies [e.g Ge, 1997; Méheust and Schmittbuhl, 2001, 2003; Boutt et al., 2006; Cardenas et al., 2009; Wang and Cardenas, 2014, 2017; Yang et al., 2019; Yoon and Kang, 2021], as well as laboratory experiments [e.g Plouraboué et al., 2000; Méheust and Schmittbuhl, 2000; Detwiler et al., 2000; Boschan et al., 2007, 2008; Ishibashi et al., 2015].

More recently, flow channeling in fractured media has been recognized as a critical control on heat transport as well [e.g Geiger and Emmanuel, 2010; Luo et al., 2016]. Based on numerical simulations of flow and heat transport, Neuville et al. [2010b] found that the heat exchange between the rock and the fluid is either enhanced or decreased in rough fractures compared to smooth fractures with equivalent mechanical apertures, depending on the fracture's morphology and aspect ratio. They concluded that because of the presence of larger flow velocities, leading to reduced transit times in the channeled areas, the heat exchange is generally less efficient in fractures with variable apertures compared to smooth fractures (the so-called parallel plate) with the same hydraulic aperture. The authors applied their modeling approach to the geothermal reservoir of Soultz-sous-Forêts, France, leading to predictions of decreased thermal exchanges in rough fractures compared to smooth ones having identical hydraulic transmissivities [Neuville et al., 2010a]. In their modelling studies, Neuville et al. [2010b,a, 2011] neglected temporal and spatial temperature variations within the rock matrix.

Neuville et al. [2013] moved beyond the assumption of constant matrix temperatures and demonstrated that the hydrothermal behavior within a fracture is heavily influenced by fracture-matrix heat exchange processes. In their study, the Navier-Stokes equations and advection-diffusion equations were solved in a simple 3-D model of a fracture consisting of flat parallel walls perturbed by a single sharp asperity [Neuville et al., 2013]. More recently, using a 3-D numerical model in which the flow in the fracture is described by the Reynolds equation, Fox et al. [2015] offered practical understanding of the effects that fracture aperture variations have on heat production in geothermal reservoirs. Notably, Fox et al. [2015] demonstrated that the thermal exchanges between the fluid flowing within a rough-walled fracture and the surrounding matrix are reduced in comparison to planar
surfaces, because flow channeling reduces the contact area over which heat conductive transfer takes place. They conclude that, as a consequence, fracture aperture variations generally have a negative impact on the thermal performance of geothermal reservoirs. More recently, from theoretical and numerical analyses of heat transfer in geometrically simple 3-D conceptual models of fractures, Klepikova et al. [2016] demonstrated that flow channeling locally enhances heat diffusion rates because a channel (cylindrical conduit) is more efficient in exchanging heat between a fracture and the matrix than a planar fracture of equivalent surface. Regardless of the different underlying assumptions, these results reveal that two opposing effects related to flow channeling impact fracture-matrix heat exchange at the fracture scale. On the one hand, heat transfer is locally enhanced by increasing the dimensionality of the diffusive flux [Klepikova et al., 2016]; on the other hand, heat transfer is reduced by decreasing the effective contact area [Fox et al., 2015; Guo et al., 2016]. Consequently, it is necessary to jointly quantify the underlying mechanisms using high-fidelity physics and modelling, in order to understand which of those opposing effects is dominant, and under which conditions.

The results of recent field experiments also call for the development of numerical modelling that considers flow channeling within individual fractures. Heat tracer experiments have been conducted recently, for example, at the experimental site of Ploemeur, France (H+ observatory network) [Read et al., 2013; Klepikova et al., 2016; de La Bernardie et al., 2019], in a field site in Altona, NY [Hawkins et al., 2017] and at the Grimsel Test Site (GTS), Switzerland [Doetsch et al., 2018]. These in situ experiments have demonstrated that, due to the signature of fracture heterogeneity on heat transfer processes, predictions based on the classical parallel plate fracture model differ significantly from field observations in terms of the first arrival time, the maximum amplitude and the tailing of the thermal breakthrough.

We present a numerical study in which the fracture is described in two dimensions (2-D) and the impermeable rock matrix in three dimensions (3-D). The flow in the fracture is described by the Reynolds equation, that is, assuming that the lubrication approximation (and hence, the local cubic law) is valid, while heat transport is described by the advection-diffusion equation in 2-D in the fracture plane, and in 3-D in the matrix. This formulation allows us to investigate heat transport along the fracture plane and in the matrix, as well as heat exchange between them, while allowing for much faster numerical simulations than with a 3-D discretization of the fracture.
We simulate 20 different rough topographies with a Hurst exponent $\zeta = 0.8$, with aperture closures $\gamma/d_m$ varying from 0.001 to 0.6 over a wide range of mean fracture apertures, and with four different values of the mismatch scale $L/L_c = 1, 2, 4, 16$. More precisely, we investigate how these properties impact heat exchange at the fluid/rock interface and heat transport along a fracture, in terms of the mean behavior among 20 fractures having such a geometry. We address how flow heterogeneity within the fracture affects the macroscopic properties (i.e., the hydraulic transmissivity, the velocity of the thermal front, and its width) ultimately governing the efficiency of the fluid mass and heat transport through the fracture.

Compared to previous studies considering simplified fracture geometry [Gringarten et al., 1975; Kocabas, 2005; Pruess and Doughty, 2010; Jung and Pruess, 2012; Neuville et al., 2013; Klepikova et al., 2016; Zhou et al., 2017] and/or a simplified heat transfer model [Neuville et al., 2010b,a, 2011], the developed numerical model of flow and heat transport considers simultaneously heat conduction in the matrix and in the fracture, as well as heat advection coupled to flow channeling in the fracture. Transient alteration of the fracture's geometry due to thermo-hydro-mechanical-chemical (THMC) coupling [Tsang, 1991; Taron and Elsworth, 2009; Pandey et al., 2014; Guo et al., 2016; Salimzadeh and Nick, 2019; Patterson and Driesner, 2020] are not considered as such effects are mainly relevant for very sharp temperature contrasts and typically act over time scales of months or years. Compared to recent works investigating heat transfer within rock samples [Luo et al., 2017; Chen and Zhao, 2020], our modelling results allow evaluating the impact of the (statistical) geometrical properties of a single geological fracture on heat transfer and to determine the key controlling parameters. Ultimately, this work aims at providing improved parameterizations and guidance for effective low-dimensional fracture models. The presented results also advance our understanding of how fracture heterogeneity control the efficiency of diffusive exchange processes at the fracture scale.

This paper is organized as follows. Section 2 describes the physical conceptualization and the implemented numerical approach. The numerical results are presented in Section 3, where we first describe the results for a given aperture-field realizations, and then present the general trends that are observed when considering large sets of synthetic fracture fields. The relevance of our models to practical configurations and applications is discussed in Section 4. Finally, Section 5 concludes with a summary of the most important findings, and outlines possible future developments.
## 2. Methods

To investigate the sensitivity of the hydrothermal behavior of a rough-walled fracture embedded in a homogeneous rock matrix to statistical fracture aperture properties, we develop a numerical model of flow and heat transport. We first present the self-affine aperture fracture model. Then, we develop a numerical model of flow and heat transport using the finite element-based software COMSOL Multiphysics® [COMSOL Multiphysics User's Guide, Version 5.4., 2018], as detailed in the following.

### 2.1. Roughness of fracture aperture

We consider fractures whose projection on the mean fracture plane is square and of lateral length $L$. Experimental studies carried out on cores from natural joints [Brown et al., 1986] have shown that the two fracture walls are self-affine, but are matched, that is they display identical topography fluctuations, at scales larger than a critical length scale $L_c \le L$. This scale, also denoted mismatch scale, is the upper limit to the self-affinity of the aperture field. It is the only characteristic scale smaller than $L$ available to describe the aperture geometry, and is a property related to the regimes of faulting (e.g., strike-slip, normal) and the history of the fracture (erosion, dissolution, precipitation processes) which is independent of the fracture size $L$ [de Dreuzy et al., 2012].

Using an algorithm adapted from Méheust and Schmittbuhl [2003], we generate fracture aperture fields $d_{fr}(x, y)$ with periodic boundary conditions that are self-affine up to $L_c$ and have a mechanical aperture $d_m$ (i.e., the distance between the mean planes of the fracture walls, which are parallel to each other; if the walls are nowhere in contact, then the mechanical aperture is the average value of $d_{fr}(x, y)$), and standard deviation $\gamma$ of the aperture field over the entire fracture. In our model, the Hurst (or, roughness) exponent has been chosen constant at $\zeta = 0.8$, a value observed for many natural and artificial fracture surfaces [Schmittbuhl et al., 1993; Bouchaud, 1997; Renard et al., 2013]. Each fracture is characterized by the fracture closure $\gamma/d_m$, which expresses the vertical extent of roughness relative to the fracture wall separation. To keep a simple boundary geometry of the domain, we prevent contact between fracture surfaces by only considering relatively small fracture closures. Consequently, the mechanical aperture is also the mean aperture. Due to stochastic variations, the mean aperture of the realizations deviates slightly from the value specified when generating the field. Here and below,
d_m refers to the actual mechanical apertures of each fracture realization, and the size of the fracture is fixed to L = 10 m. The grid size is 1024 × 1024. Using different seeds of the random generator of the white noise used in the algorithm, it is possible to generate multiple independent self-affine aperture realizations with the same underlying statistical parameters. An example of a fracture aperture profile is shown in Figure 1b.

## 2.2. Hydrothermal modelling

### 2.2.1. Hydraulic flow

Flow within a fracture is modeled as a steady-state flow where viscous forces dominate inertial effects, that is, at a low Reynolds number. Furthermore, we apply the lubrication approximation, according to which fracture walls have small local slopes [Zimmerman and Yeo, 2013]. Under these assumptions, out-of-plane (i.e., along z) components of fluid flow become negligible, and the velocity field is dominated by in-plane components. Consequently, the hydraulic flow through the fracture, $q$, defined as the integral over the local fracture aperture of the fluid velocity $u$, can be related to local apertures through a local cubic law similar to the law used for smooth (parallel plate) fractures [Méheust and Schmittbuhl, 2001; Zimmerman and Yeo, 2013]:

$$q = -\frac{d_{fr}^3(x, y)}{12\eta} \nabla P, \qquad (1)$$

where $P$ is the local pressure [Pa] and $\eta$ is the fluid's dynamic viscosity [Pa s]. $P$ and $q$ only depend on the two spatial coordinates that define the fracture's mean plane, which will be simply denoted "fracture plane" in the following. For quasi-parallel flows, the Reynolds number is generally defined as $Re = U_{charact}(\rho_f l_z^2)/(\eta l_h)$ [Méheust and Schmittbuhl, 2001; Neuville et al., 2011], where $\rho_f$ is the fluid density [kg/m³], $l_z$ and $l_h$ denote estimates of the vertical and horizontal scales of variation of the velocities [m], $U_{charact}$ is a characteristic velocity which we choose equal to the maximum velocity within a parallel plate fracture geometry of aperture $d_m$, that is, $u_M = \nabla P d_m^2/8\eta$ as estimated from the classical cubic law. In the particular case of a rough fracture, one can consider $l_z = d_m$ and $l_h = L_c$, so $Re$ can be expressed as

$$Re = u_M \frac{\rho_f d_m^2}{\eta L_c}, \qquad (2)$$
For the flow to remain linear in the fracture, the Reynolds numbers should be low, that is, $Re < 1$ [Oron and Berkowitz, 1998; Brush and Thomson, 2003; Lee et al., 2015].

Furthermore, we assume the fluid to be incompressible ($\nabla \cdot \mathbf{u} = 0$), which implies that $\mathbf{q}$ is also conservative: $\nabla \cdot \mathbf{q} = 0$ [Plouraboué et al., 2000]. Inserting the local cubic law (Equation (1)) in this conservation law yields the Reynolds equation:

$$ \nabla \cdot (d_{\text{fr}}^3(x, y)\nabla P) = 0. \qquad (3) $$

As boundary conditions, we impose the pressure at the inlet ($x = 0$) and outlet ($x = L$) of the fracture, resulting in a macroscopic pressure gradient $\nabla P$, and consider periodic boundary conditions at $y = 0$ and $y = L$. The aperture field of the fracture also has periodic boundary conditions by construction (i.e., apertures at $x = L$ are identical to those at $x = 0$). Although this condition is artificial, it is more appropriate than assigning no-flow boundaries, which could restrict the flow, thus limiting the sensitivity to fracture properties [e.g. Oding, 1992; Oron and Berkowitz, 1998; Méheust and Schmittbuhl, 2001]. The surrounding rock matrix is assumed to be impermeable. An example of the hydraulic flow computed inside a fracture's aperture field, a profile of which is shown in Figure 1b, is shown in Figure 1a as black arrows.
<figure><img src="image_4.png" /><figcaption>Figure 1: (a) 3-D sketch of a fracture model. The x-axis is along the mean hydraulic flow, the y-axis is along the mean fracture plane and perpendicular to the x-axis, and the z-axis denotes the out-of-plane direction (with respect to the mean plane). Filled color represents the temperature anomaly after 1000 s of injection. The arrows represent the direction of the flow within the fracture with the arrow length proportional to the computed flow velocity. The fracture is embedded in a homogeneous and impermeable rock matrix. The horizontal extent of the fracture is 10 m×10 m. (b) The local fracture aperture *d* profile along *y* = 5 m (shown by yellow dashed line in (a)). The mean fracture aperture *d*<sub>m</sub> is shown by the dashed line.</figcaption></figure>

The permeability of a rough fracture depends both on the mean aperture
and the geometry of the rock walls. The hydraulic aperture for a rough frac-
ture is classically defined as the aperture of the parallel plate (i.e., smooth) fracture of identical transmissivity. It can thus be computed from the measured flux $Q$ and imposed pressure gradient $\Delta P$ along the fracture [Tsang, 1992] as:

In the field, the hydraulic aperture can be assessed based on hydraulic transmissivity measurement obtained, for example, by pumping tests, while mechanical aperture can be assessed from roughness analysis on core material.

### 2.2.2. Heat transport

The 3-D finite element modeling tool COMSOL Multiphysics® allows for straightforward coupling of fluid flow and heat transport. We consider conductive and advective heat transport in the fracture, and heat conduction in the surrounding rock matrix. Moreover, our modelling approach allows avoiding 3-D discretization of a thin fracture domain, which would lead to very large computational times when considering realistic ratios of the size of the matrix along $z$ to the mean fracture aperture. Figure 1a presents a 3-D sketch of our model with a fracture located within the $x-y$ plane that is surrounded by the impermeable rock matrix at $|z| > 0$.

We assume that fluid at a constant pressure and temperature of $T_{inj}$ enters the fracture, initially at temperature $T_{rock}$, at the left model boundary ($x=0$), and flows in response to the imposed pressure gradient from left to right. The rock temperature at the outer boundaries as well as the temperature at the fracture outlet are $T_{rock}$. The injected fluid temperature is here warmer ($T_{inj} > T_{rock}$) than the initial rock temperature, a scenario typically encountered during hydrological investigations, in hydro-carbon recovery or in nuclear waste leakage. In the following, we shall consider the relative temperature deviation from the host rock's initial temperature (defined in Eq. 10), so the results for the injection of a colder fluid ($T_{inj} < T_{rock}$), a scenario typical of geothermal systems, would be exactly identical.

Several studies have found that heat conduction in the matrix in the direction parallel to the fracture has only a minor effect on the temperature distribution [e.g. Jung and Pruess, 2012]. We do consider horizontal conductive heat transport in the matrix, but assume that the conductive heat flux through the boundaries at $x, y < 0$ and at $x, y > L$ can be neglected,

$$d_h = \left( \frac{12\eta Q}{\Delta P} \right)^{\frac{1}{3}}. \qquad (4)$$
and hence we do not extend the rock matrix in these directions beyond the fracture length. The thickness of the rock matrix layer $D = 2$ m was chosen sufficiently large not to influence the temperature field within the fracture during the simulation time. The boundary conditions for temperature at the lateral boundaries of the fracture are periodic, similarly to those for the aperture field and local flux field. We neglect natural convection by temperature-induced bouyancy effects, that is, we consider that the fluid's density is independent of its temperature. The heat transport in the system can then be described as follows.

In the impermeable rock matrix:

$$ \rho_r C_{p,r} \frac{\partial T}{\partial t} + \nabla \cdot \mathbf{q}_r^{(h)} = Q^{(h)}, \qquad (5) $$

where the conductive heat flux is given by

$$ \mathbf{q}_r^{(h)} = -k_r \nabla T. \qquad (6) $$

In the fracture:

$$ d_{fr} \rho_f C_{p,f} \frac{\partial T}{\partial t} + d_{fr} \rho_f C_{p,f} \mathbf{u} \cdot \nabla_t T + \nabla_t \cdot \mathbf{q}_{fr}^{(h)} = Q^{(h)}, \qquad (7) $$

where the conductive heat flux is given by

$$ \mathbf{q}_{fr}^{(h)} = -d_{fr} k_f \nabla_t T. \qquad (8) $$

Here $\nabla_t$ denotes the gradient operator restricted to the fracture's tangential plane, $\rho$ the density [kg/m³], $C_p$ the heat capacity [J/kgK], $k$ the thermal conductivity [W/mK], $\mathbf{u}$ the local fluid velocity field [m/s], and the f, fr and r subscripts denote the fluid, fracture and rock, respectively.

The transport is characterized by the dimensionless thermal Péclet number $Pe$, which is the ratio between the characteristic times of heat diffusion/conduction and advection in the fracture [e.g. Ge, 1998; Gossler et al., 2019]:

$$ Pe = \frac{u_M d_m \rho_f C_{p,f}}{k_r}. \qquad (9) $$

The numerical model relies on a 2-D discretization of the fracture plane and the 3-D discretization of the matrix domain. To capture with sufficient accuracy the relative variations of temperature, we imposed a finer mesh size
around the fracture (~0.02 m), while the coarser elements (~0.2 m) are located near the outer boundaries of the rock matrix. We verified that refining the element size by a factor of 2 did not influence the resulting temperature field significantly (<0.1 %). Since we ignore thermo-hydro-mechanical-chemical (THMC) processes in this study (see Introduction and Discussion), all physical properties are assumed to be time-invariant.

For all the computations done in this study, the pressure gradient was chosen such that the Péclet number is Pe = 51. Such a high Péclet number implies that heat advection in the fracture is much faster than heat conduction in the matrix. The simulations are done at low Reynolds numbers, the maximum Reynolds number being Re = 0.8 (for fracture apertures as high as $d_m$=23 mm), which is compatible with the lubrication approximation (and hence, the local cubic law). Note, that while the definition of Péclet and Reynolds numbers varies between studies, the range of velocities and apertures studied herein is similar to those reported in previous works [e.g. Neuville et al., 2013, 2011]. Here and below, thermal parameters are selected to represent granite, the host formations of most deep geothermal projects (Table 1).

<table><caption>Table 1: Chosen rock properties representative of granite. The thermal parameters are taken from Incropera and DeWitt [1996]; Klepikova et al. [2016] and Kant et al. [2017].</caption><thead><tr><td>Parameter</td><td>Matrix, r</td><td>Water, f</td></tr></thead><tbody><tr><td>Density ρ, kg/m³</td><td>2500</td><td>1000</td></tr><tr><td>Thermal conductivity k, W/(m K)</td><td>0.59</td><td>3.5</td></tr><tr><td>Heat capacity C<sub>p</sub>, J/(kg K)</td><td>750</td><td>4189</td></tr></tbody></table>

## 3. Results

### 3.1. Hydraulic behaviour

In Figure 2a we present the ratio of hydraulic to mechanical apertures $d_h/d_m$, as a function of the fracture closure $\gamma/d_m$ for 20 families of fractures. A family of fractures refers here to a set of fractures generated with the same random seed, but with different fracture closure $\gamma/d_m$ and/or values for the mismatch scale $L_c$ (investigated in Section 3.2.5). In Figure 2a each curve represents a family of fractures with the same rock walls, but their separation ($d_m$) differs. The behavior of a parallel plate model corresponds to the horizontal dashed line $d_h/d_m = 1$. For closure $\gamma/d_m < 0.1$, the hydraulic
341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

behaviour is close to the parallel plate model. For different fracture families, fluid flow tends to be channelized and different hydraulic behaviors can be observed for similar values of the fracture closure $\gamma/d_m$. In agreement with previous studies [e.g. Méheust and Schmittbuhl, 2001, 2003], we find that the hydraulic behavior of rough fractures deviate monotonically from the ideal parallel plate model as fracture closure is increased. Depending on the geometry as determined by the random seed, the deviation from the parallel plate model can be positive, which implies that the fracture is more conducive to flow than a parallel plate with identical mean separation (flow-enhancing behavior) [Méheust and Schmittbuhl, 2000]. If the deviation is negative, then the fracture is characterized by flow-inhibiting behavior. Similar to previous studies [Méheust and Schmittbuhl, 2001; Neuville et al., 2010b], we see that, for most cases (75% of fracture families), the effective hydraulic transmissivity at the scale of the fracture is reduced.

In order to better understand the origin of these differences in hydraulic behaviour, we plot in Figures 2b-d some of the investigated aperture fields $d_{fr}(x, y)/d_m$ for the highest closure considered, and in Figures 2e-g the corresponding maps of local fluxes (2-D velocities) normalized by their mean value ($q/||q||$). In Figure 2b, we see a large channel oriented parallel to the applied pressure gradient (from left to right), which constitutes a configuration favorable to flow as seen in Figure 2e. This fracture morphology corresponds to the largest ratio $d_h/d_m$ in Figure 2a. Figure 2c shows the map of the ratio $d_{fr}(x, y)/d_m$ for one of the fracture families considered in Figure 2a; this family demonstrates a moderate flow-inhibiting behavior (family A, red markers). In this case, the fracture aperture field is characterized by a main tortuous channel with smaller flow obstacles (Figure 2f). In Figure 2d a barrier is seen across the whole fracture that is perpendicular to the applied pressure gradient. As shown in Figure 2g, this results in a strong flow-inhibiting behaviour of the fracture (lowest ratio $d_h/d_m$ in Figure 2a). In general, the resulting local fluxes in fractures vary over several orders of magnitude with a ratio of the local flux to the mean local flux, $q/||q||$, reaching ~12 in some geometrical configurations.

Additional simulations with other mismatch scales, $L_c$, indicate, in agreement with Méheust and Schmittbuhl [2003], that the mismatch scale also has a critical impact on the flow channeling. We find that the mean hydraulic behavior of rough-walled fractures generally converges to the parallel plate estimate when the ratio $L/L_c$ increases. These results are discussed later in relation to the resulting thermal behaviour (Section 3.2.5).

373

374

375

376

377

378
<figure><img src="image_5.png" /></figure>

Figure 2: (a) Evolution of the ratio of the hydraulic aperture to the mechanical aperture, $d_h/d_m$, as a function of the fracture closure $\gamma/d_m$ for 20 fracture families representing different fracture aperture topographies. The Hurst component is $\zeta = 0.8$ and the ratio of the fracture length to the correlation length is $L/L_c = 1$. (b-d) Maps of $d_{fr}(x, y)/d_m$; for the fracture demonstrating the largest ratio $d_h/d_m = 1.56$ at the largest closure considered, $\gamma/d_m = 0.56$ (b); for fracture family A (red markers in (a)), which exhibits a moderate flow-inhibiting behavior $d_h/d_m = 0.59$ at the largest closure considered, $\gamma/d_m = 0.6$ (c); for the fracture demonstrating the smallest ratio $d_h/d_m = 0.18$ at the largest closure considered, $\gamma/d_m = 0.56$ (d). (e-g) Maps of local fluxes (2-D velocities) normalized by the mean local flux from left to right, corresponding respectively to geometries (b-d).
3.2. Thermal behaviour

3.2.1. Thermal front definition

Figure 3 presents snapshots of the temperature field simulated using a parallel plate fracture model and a rough fracture model. Here and throughout the paper (see Table 2), we consider as leading example the self-affine fracture from family A shown on Figure 2c at closure $γ/d_m= 0.46$. Family A is considered as a representative example since this fracture family demonstrates moderate flow inhibiting behaviour as most natural rock fractures [Méheust and Schmittbuhl, 2001]. We calculate the non-dimensional temperature anomaly as

$$ \Delta T(x, y, z, t) = \frac{T(x, y, z, t) - T_{\text{rock}}}{T_{\text{inj}} - T_{\text{rock}}} \quad (10) $$

The temperature distribution in the rough-walled fracture is highly heterogeneous (Figure 3b) and the temperature evolution over time may differ considerably even for points located close to each other (Figure 3e). Initially, the thermal anomaly propagates along preferential large aperture channels and reaches for instance points B, D and C of Figure 3b, whose temperature evolution is shown in Figure 3c. At these points, the rate of change in temperature slows down after the first few tenths of seconds. A similar trend, albeit less pronounced, is observed for a flat fracture (point A in Figure 3a). On the contrary, the temperature at the points in regions of the fracture of low local aperture has a slower dynamic (point E). Overall, the variation of the temperature field over time and space is complex. As shown in Figure 3 (c) and (d), thermal plumes advance approximately 0.1 m into the rock matrix. These observations confirm that considering the thickness of the rock matrix layer $D = 2$ m is sufficient to eliminate the boundary effect.

To evaluate how the temperature field is linked to the pressure gradient, we use here the concept of a thermal front. The thermal front's velocity is an essential parameter for a geothermal reservoir as the cold front arrival associated with the (re-)injection of fluids causes a decrease of the temperature of the produced fluid, thus determining the longevity and economic prospect of the system [e.g. Nottebohm et al., 2012]. Furthermore, a widely used concept to characterize hydraulic properties from heat tracer tests is based on measuring thermal velocities, that are generally derived from thermal breakthrough curves using predefined values between injection and initial temperatures [Gossler et al., 2019]. We define the thermal front as the set of locations at which $\Delta T = 1/2$ (black lines in Figure 3).
414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

For both the parallel plate and heterogeneous fractures considered above, we find as expected that heat loss at the fracture walls creates a thermal front (black solid line in Figure 3) that is delayed relative to the fluid front (blue solid line in Figure 3) [e.g. Bodvarsson, 1972]. Thus, for the parallel plate fracture, when the fluid front is approaching the outlet of the fracture, the thermal front travels a normalized distance $x_{PP}^*/L$ equal to 0.3 (Figure 3a).

In order to characterize how the thermal behavior evolves on average, we consider the evolution of the thermal front with time. Figure 3b presents an example of the thermal front advancement for the fracture family A, $\gamma/d_m = 0.46$ (black dashed lines). The front spreading pathway varies with the local fluid velocity due to the roughness of the fracture aperture. As the thermal front grows, heat fingers are developing along preferential flow paths within the fracture plane, mainly in the middle region of the fracture. This causes deviations of the thermal front position from its average $x_{ROUGH}^*/L$. In the following, we characterize the advancement and the evolution of thermal fronts in rough fractures and determine geometrical parameters of individual fractures controlling this advancement. This is achieved by studying the hydrothermal behavior of models with different fracture aperture patterns. Furthermore, we compare the results with the reference case of a fracture modeled with two parallel plates separated by a constant aperture $d_m$ (i.e., no self-affine spatial variations). The key geometric characteristic of the studied fractures are presented in Table 2.
<figure><img src="image_6.png" /></figure>

Figure 3: Comparison of temperature anomalies in a parallel plate fracture ((a) - fracture plane view and (c) - view of the mid-longitudinal cross section) and in a model with strong fracture roughness belonging to fracture family A (see Figure 2c), with $\gamma/d_m = 0.46$ and $L/L_c = 1$ ((b) - fracture plane view and (d) - mid-longitudinal cross-section view). The pressure gradient was chosen such that the Péclet number was the same for both simulations, $Pe = 51$. Both fractures have the same mechanical aperture. The fluid and thermal front profiles, shown by blue and black solid lines, respectively, are obtained at $t = 900$ s, that is slightly below the time $t_{sh} = 920$ s necessary for a volume equal to the total volume of the fracture to flow through the fracture. Dashed black lines depict three thermal front profiles obtained at different times ($t = 10, 100, 500$ s). The letters A-E indicate the locations where the temperature evolution is observed in Figure (e). (e) Temperature as a function of time at location A for a parallel plate fracture (red line) and at locations B-E for the heterogenous fracture (black lines).
<table><caption>Table 2: Geometric parameters of fractures studied in this work.</caption><thead><tr><th>Test ID</th><th>Fracture family</th><th>Closure</th><th>Mechanical aperture</th><th>Correlation</th></tr><tr><td></td><td></td><td>γ/d<sub>m</sub></td><td>d<sub>m</sub>, mm</td><td>length L/L<sub>c</sub></td></tr></thead><tbody><tr><td>Test 1</td><td>family A</td><td>0.001–0.6</td><td>2.5 × 10<sup>-3</sup>–1.5</td><td>1</td></tr><tr><td>Test 2</td><td>equivalent par-allel plate</td><td>–</td><td>2.5 × 10<sup>-3</sup>–1.5</td><td>–</td></tr><tr><td>Test 3</td><td>family A</td><td>0.05–0.6</td><td>1.25–15</td><td>1</td></tr><tr><td>Test 4</td><td>equivalent par-allel plate</td><td>–</td><td>1.25–15</td><td>–</td></tr><tr><td>Test 5</td><td>family A</td><td>0.32–0.6</td><td>12.5–23</td><td>1</td></tr><tr><td>Test 6</td><td>equivalent par-allel plate</td><td>–</td><td>12.5–23</td><td>–</td></tr><tr><td>Test 7</td><td>20 families</td><td>0.05–0.6</td><td>1.25 × 10<sup>-1</sup>–1.5</td><td>1</td></tr><tr><td>Test 8</td><td>family A</td><td>0.02–0.6</td><td>5 × 10<sup>-2</sup>–1.5</td><td>2, 4, 16</td></tr></tbody></table>

### 3.2.2. Influence of fracture closure

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

We now investigate how the roughness amplitude influences the thermal behaviour of fractures from family A (Figure 2c), which exhibits a moderate flow-inhibiting behavior (red markers in Figure 2a). To do so, we vary the fracture closure $γ/d_m$ by varying the fracture’s mechanical aperture $d_m$, while keeping the same standard deviation for their height distributions, $γ$ (Test 1, Table 2). Examples of fracture apertures generated on a 1024 × 1024 grid are shown in Figure 4a. For small fracture closure $γ/d_m = 0.02$ (Figure 4a top), spatial variations of the aperture are negligible in comparison to the mean aperture. As the fracture is closed, $γ/d_m = 0.21, 0.40$ and $0.59$ (Figure 4a from top to bottom), relative fluctuations of the aperture increase.

As a consequence of increased closure, flow channeling becomes more important. This is shown in Figure 4b, where maps of local fluxes normalized by the mean local flux are shown. In Figure 4b, for high closure cases, the flow tends to avoid regions of small local apertures and, consequently, is localized in a large aperture channel along the flow direction. This channel can be seen across almost the whole fracture in Figure 4a (bottom maps). The aperture of this channel is relatively small in the vicinity of the inlet and in the vicinity of the outlet, leading to the flow-inhibiting behavior displayed in Figure 2. Finally, Figure 4c presents the simulated temperature fields
<footnote n="457">457 (after 1000 s of injection) for different values of fracture closure.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 4: (a) Maps of the ratio $d_{fr}(x, y)/d_m$ for $\gamma/d_m$ = 0.02, 0.21, 0.40, and 0.59, for fracture family A. Blue colors indicate areas of small local aperture, while orange shadings denote larger apertures. (b) Maps of local fluxes (2-D velocities) normalized by the mean local flux. (c) Temperature anomaly induced after 1000 s of injection, which approximately corresponds to the time necessary to replace the total volume of the fracture. The thermal fronts are shown as black lines; the fluid flows from left to right and the hot fluid is injected from the left. The linear color scale is the same for all figures of each column.</figcaption></figure>
458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

As discussed above, one of the important characteristics of the geometry of the mixing zone (where $0 < \Delta T < 1$) is the position and shape of the thermal front (black dashed lines in Figure 4). The shape of the thermal front within the fracture is strongly correlated with the hydraulic flow. For small values of fracture closure, the thermal front is almost straight and transverse to the mean flow direction. However, the thermal front becomes less smooth as $\gamma/d_m$ increases, the pattern of temperature distribution becomes complex with 'slow zones' forming in regions of low local fracture apertures and thermal fingers developing along preferential flow channels. Thus, for $\gamma/d_m$=0.59, when the most advanced thermal finger passes half of the fracture's length, the most delayed region of the thermal front are still in the vicinity of the fracture inlet (Figure 4c bottom). Hence, the width of the thermal front parallel to the flow direction increases as the fracture is closed.

In order to quantify the variability of thermal front velocities, we use the results of Test 1 (Table 2) and observe how the mean position of the thermal front $x_{ROUGH}^*$ evolves with time. The standard deviation of the front position, which quantifies its width along the mean flow direction, is plotted in Figure 5 against the position $x_{ROUGH}^*$, at positions $y = [0, 0.01, 10]$ m. When the roughness amplitude increases, the standard deviation of the front position increases, implying that the thermal channeling effect is more pronounced, as expected. We also note that the velocity of the thermal front decreases as the fracture closure increases. The latter is related to the hydraulic behaviour of the fracture, which tends to inhibit the hydraulic flow as shown in Figure 2a.
<figure><img src="image_8.png" /><figcaption>Figure 5: Mean position of the thermal front normalized by the fracture size, $x_{\text{ROUGH}}^{*}/L$, versus the standard deviation of the front position normalized by the fracture size, $\sigma_{x\text{ROUGH}}/L$, for various fracture closure $\gamma/d_{\text{m}}$. Fracture family A, fracture closure $\gamma/d_{\text{m}}$= 0.001, 0.03, 0.10, 0.21, 0.32, 0.46, 0.60 and $L/L_{\text{c}} = 1$. Blueish marker colors refer to open fractures, and magenta marker colors refer to more closed fractures.</figcaption></figure>

We further verified this observation by comparing the results of hydrothermal simulations in rough-walled fractures with simulations in parallel plate fractures of identical mechanical aperture (Test 2, Table 2). Figure 6 presents the evolution in time of the average velocity of the thermal front $v_{\text{ROUGH}} = x_{\text{ROUGH}}^{*}/t$ relative to the $v_{\text{PP}} = x_{\text{PP}}^{*}/t$, that is, the thermal front velocity in a fracture modeled with a constant aperture $d_{\text{m}}$. For a small fracture closure $\gamma/d_{\text{m}}$= 0.05, which corresponds to a nearly smooth aperture field, a parallel plate model reproduces a similar thermal profile and the ratio $v_{\text{ROUGH}}/v_{\text{PP}}$ is close to 1. As the fracture closure is increased, $\gamma/d_{\text{m}}$ = 0.18, 0.32, 0.46
and 0.60 (Figure 6), the velocity of thermal front becomes slower compared to thermal front velocity in parallel plate fractures with identical mechanical aperture $v_{ROUGH} < v_{PP}$.

<figure><img src="image_9.png" /><figcaption>Figure 6: Time evolution of the average velocity of the thermal front normalized by the thermal front velocity in the equivalent parallel plate fracture (i.e, the parallel plate of aperture equal to the rough fracture’s mechanical aperture $d_m$), for fracture family A, fracture closures $\gamma/d_m = 0.05, 0.18, 0.32, 0.46, 0.6$, and $L/L_c = 1$. Square markers refer to fractures with large mechanical apertures $d_m$, and circle markers refer to smaller apertures $d_m$.</figcaption></figure>

For all simulations, higher ratios of thermal front velocities $v_{ROUGH}/v_{PP}$ are observed close to the fracture inlet. Figure 6 also demonstrates that for small fracture closure $\gamma/d_m = 0.05$, $v_{ROUGH} > v_{PP}$, meaning that the thermal
front advances faster in rough fractures compared to what would be expected with a parallel plate fracture of (uniform) aperture $d_m$. However, after a short transient regime, the thermal front velocities in rough and in smooth fractures become equal, $v_{ROUGH} = v_{PP}$. As the thermal tracer enters the fracture, heat diffuses first across the fracture, and, afterwards, away into the matrix. Once the rock temperature starts to evolve in time and in space, the ratio $v_{ROUGH}/v_{PP}$ stabilizes. Larger fracture apertures imply that heat needs more time to diffuse across the fracture aperture. To verify this, we evaluate through Tests 3-4 (Table 2) the thermal front velocity within the fractures in family A with different mechanical apertures but the same fracture closures $\gamma/d_m$ (Figure 6). For a fracture with an aperture $d_m = 15$ mm, the thermal front velocity ratio becomes quasi-steady after $t = 500$ s, when the thermal front has travelled a mean normalized distance $x_{ROUGH}^*/L$ equal to 0.17. For a fracture with the same closure, $\gamma/d_m = 0.05$, but smaller aperture, $d_m = 1.5$ mm, the thermal front velocity ratio becomes quasi-steady after $t = 200$ s, and for fractures with small apertures $d_m < 1$ mm, the thermal front velocity ratio stabilizes already after $t = 100$ s, when the thermal front has travelled a mean normalized distance $x_{ROUGH}^*/L$ of less than 0.1. We further evaluated through Tests 5-6 (Table 2) that for more closed fractures (asterisk in Figure 6) with large apertures ($d_m = 13 - 23$ mm), the thermal front velocity ratio becomes quasi-steady after $t = 400 - 700$ s, respectively.

### 3.2.3. Influence of the fracture hydraulic aperture

For the same cases as in Figure 6, Figure 7 presents the ratio $v_{ROUGH}/v_{PP}$ versus the ratio between the hydraulic and mechanical apertures $d_h/d_m$ for two different times: for a very short duration $t = 30$ s, when the regime is still transitory (grey markers), and for a longer duration, $t > 700$ s (black markers), when the thermal front velocity ratio becomes quasi-steady. For a short duration, we can observe that for a rough aperture, the thermal front advances systematically faster in comparison to what we expect from the hydraulic behavior (all the grey points are above the 1:1 line). Our results also demonstrate that this effect is more pronounced for fractures with large apertures (circle markers referring to small apertures are closer to the 1:1 line than square markers and asterisk referring to larger apertures). The demonstrated faster propagation of the thermal signal in rough fractures ($v_{ROUGH}$) when compared to that in smooth fractures ($v_{PP}$) agrees with the work of Neuville et al. [2010b], who attributes this effect to a decrease in heat exchange efficiency in rough fractures due to the reduction of transit times
534 in the channeled areas of a rough-walled fracture.

535 However, once the hosting rock temperature starts to evolve, a process
536 which was not accounted for in the model of Neuville et al. [2010b], the ther-
537 mal front velocity in rough fractures slows down (Figure 6), and, for t > 400
538 s, we obtain (Figure 7) a perfect correlation between the ratio v<sub>ROUGH</sub>/v<sub>PP</sub>
539 and the ratio between the hydraulic and mechanical apertures d<sub>h</sub>/d<sub>m</sub> (black
540 triangle markers). This implies that once heat has diffused along the out-
541 of-plane direction over the fracture's aperture, the thermal behavior of a
542 rough-walled fracture is determined by its effective hydraulic transmissivity.
543 For the morphology investigated here (Family A fractures), the permeability
544 is reduced, and, thus, the advancement of the thermal front is slower in rough
545 fractures compared to what would be expected with a model of flat fractures
546 having the same mechanical aperture d<sub>m</sub> (v<sub>ROUGH</sub>/v<sub>PP</sub> < 1).
<figure><img src="image_10.png" /><figcaption>Figure 7: Ratio of the average velocity of the thermal front relative to that in the equiva-
lent parallel plate model (same mechanical aperture) versus d_h/d_m, for fracture family A
(L/L_c = 1).</figcaption></figure>

We now leave the specifics of family A and consider all fracture fami-
lies in Figure 2a for different closure values (Test 7, Table 2). Statistical
thermal results presented in Figure 8 confirm the near-perfect correlation
between the ratio v<sub>ROUGH</sub>/v<sub>PP</sub> and the ratio between the hydraulic and me-
chanical apertures d<sub>h</sub>/d<sub>m</sub>. This means that, once stabilized, such that heat
has diffused along the out-of-plane direction over the fracture’s aperture, the
mean position in time of the thermal front within a rough-walled fracture
is determined by the hydraulic aperture d<sub>h</sub>. For different fracture families
with the same fracture closure γ/d<sub>m</sub>, we find that the flow inhibiting be-

547

548

549

550

551

552

553

554

555
556

557

558

559

560

561

562

563

564

havior is favored statistically [Méheust and Schmittbuhl, 2001]. Note that these flow-enhancing or flow-inhibiting behaviors of individual fractures are related to the fractures’ hydraulic anisotropy, as a flow-enhancing fracture becomes flow-inhibiting (and vice-versa) when the flow direction is rotated by 90° (see discussion in Méheust and Schmittbuhl [2001]). Hence, as the fracture, due to the roughness of its walls, is either less or more permeable than a flat parallel plate of identical mechanical aperture, the efficiency in transferring heat is also highly variable (0.1 < v<sub>ROUGH</sub>/v<sub>PP</sub> < 1.6) from one fracture family to another, but ratios smaller than 1 are favored statistically.
<figure><img src="image_11.png" /></figure>

Figure 8: Grey markers: averages of the normalized thermal front velocities $v_{ROUGH}/v_{PP}$ versus the ratio of hydraulic to mechanical apertures $d_h/d_m$, for more than 150 cases whose hydraulic apertures are presented in Figure 2 for various closure $\gamma/d_m$ values, $\zeta = 0.8$, $L/L_c = 1$. Red markers correspond to fracture family A and different mismatch scales $L/L_c = 2, 4, 16$, respectively. The normalized thermal front velocities are considered at times larger than the time needed for heat to diffuse across the fracture aperture). Markers are transparent to highlight the density of points clustered near the 1 : 1 line, which holds for parallel plates separated by $d_h$. The minimum root mean square error (RMSE) used to quantify the goodness-of-fit is 0.01.

### 3.2.4. Conductive heat flux

While previous studies characterized the heat exchange efficiency of rough fractures through the use of temperature metrics [Neuville et al., 2010a, 2013; Fox et al., 2015], in this study, we provide some insights into diffusive exchange processes at the fracture scale. Using the same data (Tests 1-2, Table
2), we now calculate the total conductive heat flux between the fracture and the rock matrix for family A for different values of the fracture closure, $\gamma/d_m = 0.05, 0.18, 0.32, 0.46, 0.60$. Our results, presented in Figure 9, demonstrate that for all cases investigated here, the conductive heat flux is greater for the equivalent parallel plate fracture (i.e., the parallel plate with an aperture equal to the rough fracture's mechanical aperture $d_m$):
$$Q_{ROUGH}^{(h)} < Q_{PP}^{(h)}.$$ Moreover, the ratio $Q_{ROUGH}^{(h)}/Q_{PP}^{(h)}$ converges in time to a plateau (Figure 9, inset). Interestingly, Figure 9, showing the plateau value versus the ratio of the hydraulic to mechanical apertures $d_h/d_m$, demonstrates that the conductive heat flux between the rough fracture and the surrounding rock can be predicted from the equivalent parallel plate model. Overall, for fracture family A, both the heat flux along the fracture and the conductive heat fluxes between the fracture and the embedding rock decrease when the roughness amplitude increases. This analysis suggests that a major cause of the observed slowing down of the thermal front in rough fractures compared to smooth fractures having the same mechanical aperture (Figure 7) is related to the flow-inhibiting behavior of fractures from family A, rather than to an increase in the efficiency of the conductive heat exchange between the fluid and the rock matrix. Finally, this result confirms that for the hydrothermal conditions studied here, when the influence of heat advection in the fracture dominates the influence of conduction in the matrix, the hydraulic aperture governs the fracture's thermal behaviour. As we shall discuss in section 4, such conditions are actually relevant for most applications (heat tracer testing and geothermal systems).
<figure><img src="image_12.png" /></figure>

Figure 9: Steady state ratio of the total conductive heat flux between the fracture and the rock matrix to the same quantity measured in the equivalent parallel plate model (same mechanical aperture as the rough fracture) versus $d_h/d_m$, for fracture family A and $L/L_c = 1$. The RMSE is equal 0.005. The inset illustrates the evolution in time of that ratio before the thermal exchange quasi-steady state is reached.

### 3.2.5. Influence of the mismatch scale/correlation length

Finally, we investigate how the mismatch scale (i.e., correlation length) influences the thermal behaviour of the fracture. To do so, we modify not only the fracture closure $\gamma/d_m$, but also the ratio $L/L_c$, while keeping the same numerical seed for the generation of the rough topographies, and the same standard deviation for their height distributions (Test 8, Table 2). Figure 10 presents the simulated temperature fields (after 1000 s of injection)
601

602

603

604

605

606

607

608

609

610

611

612

613

for different mismatch scales, $L/L_c = 2$ (Figure 10a), $L/L_c = 4$ (Figure 10b), $L/L_c = 16$ (Figure 10c), and for different values of fracture closure, $\gamma/d_m = 0.02, 0.25, 0.48, \text{ and } 0.60$ (from top to bottom). Similarly to $L/L_c = 1$ (Figure 4), the thermal channeling effect is all the more pronounced as the fracture is more closed (Figure 10). For $L/L_c = 2$ (Figure 10a), the thermal front displays large scale distortions representing a large fraction of its total width. Furthermore, we observe the refinement of filaments and global flattening of thermal fronts (i.e., decrease in their roughness amplitude) with increasing $L/L_c$ (Figure 10 from left to right), which is determined by how the decrease of the mismatch scale impacts the advecting velocity field, reducing the typical scale of its heterogeneities and the 2-D velocity contrast between preferential flow paths and regions of lower velocities.
<figure><img src="image_13.png" /><figcaption>Figure 10: Comparison of temperature anomalies induced after 1000 s of injection for different mismatch scales (a) $L/L_c = 2$, (b) $L/L_c = 4$ and (c) $L/L_c = 16$ and for different fracture closures: $\gamma/d_m = 0.02, 0.25, 0.48$, and 0.60.</figcaption></figure>
614

615

616

617

618

619

620

621

622

Despite the changes in the *q*-field with increasing *L*/*L*<sub>c</sub>, we find similarly to what is observed for *L*<sub>c</sub> = *L* (black markers in Figure 8), that the fluid flow scaling translates into that of heat transport. Indeed, red markers in Figure 8 show strong correlation between the ratio *v*<sub>ROUGH</sub>/*v*<sub>PP</sub> and the ratio between the hydraulic and mechanical apertures *d*<sub>h</sub>/*d*<sub>m</sub> for *L*/*L*<sub>c</sub> = 2, 4, and 16. This result confirms that the hydraulic behaviour allows predicting thermal channeling effects and the related thermal behavior for geological rough fractures for various values of the fracture closure and various values of the ratio *L*/*L*<sub>c</sub>.

623

# 4. Discussion

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

We find that the thermal behavior of horizontal rough-walled fractures can be predicted from their hydraulic behavior, as demonstrated for a wide range of thermal Péclet numbers 6 < *Pe* < 200. Of course, this finding is expected to be all the more valid for larger Péclet numbers (*Pe* > 200). For lower Péclet values, slight deviations from this general finding was observed for *Pe* = 5. Considering lower Péclet numbers would be very demanding in terms of computational resources, but it would also not be so relevant for applications. Indeed, Péclet numbers typically take values in the range of 10-7 000 in fractured geothermal systems under production [Horne and Rodriguez, 1983; Geiger and Emmanuel, 2010; Neuville et al., 2010a], and the Péclet numbers of thermal tracer tests range between 2 and 70 000 [Klepikova et al., 2016; Hawkins et al., 2017; de La Bernardie et al., 2019], as summarized in Table 3. This suggests a broad applicability of our inferred relationships between heat transport and hydraulic behavior. Interestingly, similar conclusions have also been previously achieved for solute transport: applying the equivalent aperture size calculated based on the equivalent permeability of the system provides an acceptable prediction of solute transport [Nick et al., 2011].
<table><caption>Table 3: Péclet numbers considered in previous studies on heat transport in fractured media.</caption><thead><tr><td>Study</td><td>Scale</td><td>Péclet number</td></tr></thead><tbody><tr><td>Horne and Rodriguez [1983]</td><td>Fracture network</td><td>10-100</td></tr><tr><td>Ge [1998]</td><td>Single fracture</td><td>5-100</td></tr><tr><td>Geiger and Emmanuel [2010]</td><td>Fracture network</td><td>14-145</td></tr><tr><td>Neuville et al. [2010a]</td><td>Single fracture</td><td>7 000</td></tr><tr><td>Neuville et al. [2013]</td><td>Single fracture</td><td>10</td></tr><tr><td>Klepikova et al. [2016]</td><td>Single fracture</td><td>5 000</td></tr><tr><td>Hawkins et al. [2017]</td><td>Single fracture</td><td>2-8</td></tr><tr><td>de La Bernardie et al. [2019]</td><td>Single fracture</td><td>70 000</td></tr></tbody></table>

In our model, fracture flow and heat transport are described in 2-D. This allows us to solve transport both in the fracture and in the rock matrix, with an extent of the matrix along the direction normal to the fracture plane that is sufficiently large to avoid boundary effects. This wouldn't be possible if we had to discretize the fracture aperture in a 3-D mesh within the fracture to account for fracture flow and transport in 3-D, or it would be so demanding on computer resources that we wouldn't be able to consider ensemble statistics of fractures with identical geometrical parameters. Due to this choice, however, 3-D flow effects cannot be accounted for in our model. A number of studies have addressed such effects, and concluded that they might be significant. These include nonlinearities in the flow induced either by local tortuosity and roughness, or by inertial effects [e.g Ge, 1997; Brush and Thomson, 2003; Wang et al., 2015, 2020; He et al., 2021]. The latter are not relevant to our study, since we consider creeping flow. However, solving the flow from the Reynolds equation (i.e. the traditional local cubic law), which cannot account for tortuosity effects in the third dimension, can result in overestimation of the transmissivity (or hydraulic aperture). Nevertheless, in laminar flow regimes, contributions of these effects do not impact flow and heat transport significantly [e.g. Brush and Thomson, 2003; Lee et al., 2015].

The fracture's dimensionality may also impact heat transport through it. We consider a horizontal fracture, and thus do not need to account for buoyancy effects (resulting from the temperature-dependence of the fluid's density), which, in a subvertical fracture, may lead to convective fluid circulation [Patterson et al., 2018]. Note that some studies solving advection-diffusion equations for heat transport in three dimensions report the emergence of
667

668

669

670

671

672

673

674

675

676

677

678

679

680

681

682

683

684

3-D effects. Thus, the presence of recirculation and stagnant zones related to highly variable morphology of the fluid-rock interface may locally, within the asperity, modify the heat exchange [Andrade et al., 2004; Neuville et al., 2013]. Another possible thermal effect was studied by Klepikova et al. [2016] and de La Bernardie et al. [2019] and emerges when large local slopes (angle between the orientation of the local fracture wall and that of the fracture plane) exist, and the heat flux in between the matrix and the fracture, occurring dominantly in the direction perpendicular to the fracture walls, is oriented locally at a large angle with respect to the fracture plane. This effect was shown to have a significant impact on the scaling of heat recovery in both space and time, and the findings were supported by field experiments performed on the fractured rock site of Ploemeur, where high aperture channels (around a few cm) participate to the transport of heat [Klepikova et al., 2016; de La Bernardie et al., 2019]. Still we do not expect that accounting for 3-D effects of fracture geometry would have a significant impact on our results. In contrast to Ploemeur field site, fractures walls in this study are assumed to have small local slopes (this is essentially what the lubrication approximation means).

685

686

687

688

689

690

691

692

693

694

695

The impact of the spatial resolution of the aperture roughness has also been investigated. Additional simulations reported in the Appendix demonstrate that using slightly lower spatial resolution (i.e., downsampling the aperture field by a factor up to 8) does not significantly modify the results. These results are in general agreement with the studies of Méheust and Schmittbuhl [2001] and Neuville et al. [2011]. The former demonstrated that the Fourier modes of the aperture field corresponding to the largest scales control flow channeling in the fracture plane for the most part, and therefore, the fracture's hydraulic aperture, while the latter confirmed this finding and further demonstrated that it also holds for heat transport in rough fractures.

696

697

698

699

700

701

702

703

704

We do not consider the thermal stress acting on preferential flow paths, which reduces the effective compressive stress along these paths and, thereby, further exacerbates flow channeling, thus impacting heat exchange processes [Guo et al., 2016; Salimzadeh et al., 2018; Patterson and Driesner, 2020]. As demonstrated by recent works of Guo et al. [2016]; Patterson and Driesner [2020], the magnitude of the changes due to thermo-mechanical effects is comparable with the magnitude of the initial aperture variation. We have chosen to ignore such thermal–hydraulic–mechanical–chemical (THMC) couplings for two reasons. First, we sought to quantify the impact of fracture
surface roughness on heat transport and exchange with the rock matrix, and wanted to discriminate these purely geometrical effects from those related to thermo-mechanical couplings. Second, since THMC processes (i) arise in Enhanced/Engineered Geothermal Systems (EGS) when strong contrasts in temperature exist and (ii) are relatively slow (effects on reservoir performance are noticeable over a time scale of months-years) [Pandey et al., 2014; Guo et al., 2016; Salimzadeh et al., 2018]. Consequently, such effects are expected to be negligible within the typical duration of heat tracer tests, which is our prime target application. Introducing THMC couplings in the model is expected to produce stronger channeling and consequently a larger deviation of the fracture's hydraulic behavior from that of the smooth fracture (parallel plate) of aperture equal to the rough fracture's mean aperture. Since this would result in a change in fracture surface topography, we would still expect the main finding of the present study to hold, namely, that heat transport behavior can be predicted from the hydraulic behavior of the fracture.

## 5. Conclusions

We have investigated numerically the influence of the statistical properties of the aperture field and upscaled hydraulic behavior on heat transport in rough rock fractures with realistic geometries. The flow regime was assumed to be laminar and at low Reynolds number, and the gradient of the aperture field to be small (lubrication approximation), so that flow in the fracture could be modeled by the Reynolds equation in the 2-D fracture plane. We considered a regime where heat transport in the fracture is moderately dominant with respect to heat conduction in the rock matrix ($Pe = 51$), a configuration which is relevant for practical situations at geothermal sites and for forced hydraulic conditions usually adopted during field heat tracer tests. We analyzed 20 rough topographies with a Hurst exponent $\zeta = 0.8$, with aperture closures $\gamma/d_m$ varying from 0.001 to 0.6 over a wide range of mean fracture apertures, and with four different values of the mismatch scale $L/L_c = 1, 2, 4, 16$.

At fixed fracture closure, the deviation from the parallel plate model increases as $L_c$ is decreased. When closing the fracture, the deviation of the hydraulic and thermal behaviors from the equivalent parallel plate model increases. In general, the thermal behaviour is highly variable among a population of fractures with identical geometrical parameters. In comparison to a fracture of uniform aperture equal to the rough fracture's mechanical
aperture, 75% of the considered fracture aperture fields exhibit a slower dis- 741 placement of the thermal front along the fracture and less thermal exchange 742 between the fracture and the surrounding rock. 743

Our main finding is that under the considered conditions, thermal be-
haviour of rough-walled rock fractures only depends on the hydraulic prop-
erties. A similar conclusion was reached by Neuville et al. [2010b], who found
that the hydraulic aperture is a better parameter than closure to assess the
thermal exchange efficiency of rough fractures. In stark contrast to Neuville
et al. [2010b], we found that the heat transport along rough-walled fractures
was similarly efficient as a parallel plate (i.e., smooth) fracture of identical
hydraulic transmissivity. The thermal fronts in rough fractures are initially
slightly more advanced (at identical times) than thermal fronts in flat frac-
tures with equivalent permeabilities, but this holds only for a very short time
(i.e., the time needed for heat to diffuse along the out-of-plane direction over
the fracture's aperture). Depending on the mean fracture aperture, this tran-
sition period lasts for tens to a few hundreds of seconds, during which the
thermal front travels a mean normalized distance $x_{ROUGH}^*/L$ equal to $\sim 0.1$.
By accounting for fracture-matrix heat exchange by transverse diffusion, a
process which was neglected by Neuville et al. [2010b], the thermal behavior
of a rough-walled fracture is found to be controlled by its hydraulic aper-
ture and boundary conditions. This striking novel finding results from an
improved description of the coupled flow and heat transport.
762

The practical implication of our finding is that thermal exchanges at the
scale of a single fracture is controlled by the effective hydraulic transmissiv-
ity. Provided that thermal properties of the host rock are known, this implies
that (1) heat tracer tests are reliable for inferring effective fracture transmis-
sivity, and (2) the geothermal efficiency can be computed at field sites using
hydraulic characterization alone. Furthermore, as long as the considered
time scale doesn't allow for significant THM(C) coupling to take place, it
follows that the temporal evolution of the geothermal efficiency can be pre-
dicted over significantly large time scales using well-known low-dimensional
hydraulic parameterizations in terms of effective hydraulic properties.
763
764
765
766
767
768
769
770
771
772

Future work could address non-Stokes flow conditions (i.e., laminar but
non-linear flow) in the fracture. Another interesting prospect is the over-
all large-scale heat transport behavior in a fractured geological formation.
It is expected to depend on the combined effects of both the local scale
heterogeneity of individual fractures and the heterogeneity at the scale of a
discrete fracture network (DFN) consisting of multiple intersecting fractures.
773
774
775
776
777
778
779 The study of coupled flow and heat transport in such DFNs will be the topic
of future work.

780
781

782

# Appendix: Impacts of the spatial resolution on the hydraulic and thermal results

783

784

785

786

787

788

789

790

791

792

793

794

795

We use the same finite-difference numerical scheme as described previ-
ously to compute the flow and temperature fields in downsampled apertures,
that is, the aperture field $d_{\text{fr}}(x, y)$ is only considered at every $n$-th point of
the grid, where $n = 2, 4, 8$. In Figure 11a, the relative difference between hy-
draulic aperture $d_h$ for the downsampled apertures and the hydraulic aperture
with full resolution of the geometric aperture, $d_{h,\text{ref}}$, is evaluated. All com-
puted hydraulic apertures are closer than 0.4 per cent to the full resolution
$d_{h,\text{ref}}$. In general, our results reveal that $d_h$ is overestimated with respect to
$d_{h,\text{ref}}$ when decreasing the spatial resolution. The relative difference between
the average velocity of the thermal front calculated inside the downsampled
apertures and apertures with full resolution is shown on Figure 11b for the
case of the largest closure considered, $\gamma/d_m = 0.6$. Here, the precision of the
approximation is better than 0.6 per cent.
<figure><img src="image_14.png" /><figcaption>Figure 11: (a) Relative errors of predicted hydraulic apertures $d_h$ when downsampling the aperture data by a ratio 2, 4, and 8; the errors are plotted as a function of the fracture closure. Fracture family A (see Figure 2c), $L/L_c = 1$. (b) Relative errors made on the average velocity of the thermal front $v_{ROUGH}$, plotted as a function of the downsampling ratio. Fracture family A, with $\gamma/d_m = 0.6$ and $L/L_c = 1$.</figcaption></figure>
796 Notation

797 Symbol list

<table><tr><td>C<sub>p</sub></td><td>heat capacity</td></tr><tr><td>d<sub>fr</sub>(x, y)</td><td>fracture aperture field</td></tr><tr><td>d<sub>m</sub></td><td>mechanical aperture</td></tr><tr><td>d<sub>h</sub></td><td>hydraulic aperture</td></tr><tr><td>D</td><td>thickness of the rock matrix</td></tr><tr><td>k</td><td>thermal conductivity</td></tr><tr><td>L</td><td>fracture length</td></tr><tr><td>L<sub>c</sub></td><td>mismatch scale (correlation length)</td></tr><tr><td>P</td><td>local pressure</td></tr><tr><td>Pe</td><td>Péclet number (Pe = u<sub>M</sub>d<sub>m</sub>ρ<sub>f</sub>C<sub>p,f</sub>/k<sub>r</sub>)</td></tr><tr><td>q</td><td>local flux or 2D velocity, i.e., the integral over the local fracture aperture of u</td></tr><tr><td>Q</td><td>flux through the fracture</td></tr><tr><td>Q<sup>(h)</sup></td><td>conductive heat flux between fracture and rock matrix, over the entire fracture walls' areas</td></tr><tr><td>Re</td><td>Reynolds number (Re = d<sub>m</sub>u<sub>M</sub>/η)</td></tr><tr><td>t</td><td>time</td></tr><tr><td>T(x, y, z, t)</td><td>temperature</td></tr><tr><td>ΔT(x, y, z, t)</td><td>non-dimensional temperature anomaly, ΔT(x, y, z, t) = T(x, y, z, t)-T<sub>rock</sub>/(T<sub>inj</sub>-T<sub>rock</sub>)</td></tr><tr><td>T<sub>inj</sub></td><td>injection temperature</td></tr><tr><td>T<sub>rock</sub></td><td>initial temperature</td></tr><tr><td>u</td><td>three-dimensional fluid velocity in the fracture</td></tr><tr><td>u<sub>M</sub></td><td>maximum velocity within a parallel flat wall fracture</td></tr><tr><td>v</td><td>average velocity of the thermal front</td></tr><tr><td>x</td><td>distance along applied pressure gradient</td></tr><tr><td>x*</td><td>mean position of the thermal front</td></tr><tr><td>y</td><td>distance normal to applied pressure gradient</td></tr><tr><td>z</td><td>out of fracture plane distance</td></tr><tr><td>γ</td><td>standard deviation of the aperture field</td></tr><tr><td>η</td><td>dynamic viscosity</td></tr><tr><td>ρ</td><td>density</td></tr><tr><td>ζ</td><td>the Hurst (or, roughness) exponent</td></tr></table>

798
799 *Subscripts*

<table><tr><td>f</td><td>fluid</td></tr><tr><td>fr</td><td>fracture</td></tr><tr><td>r</td><td>rock</td></tr><tr><td>ROUGH</td><td>rough fracture</td></tr><tr><td>PP</td><td>parallel plate fracture</td></tr></table>

800

801 *Acknowledgments*

802 Maria Klepikova acknowledges the Equal Opportunity postdoctoral funding from the Faculty of Geosciences and Environment, University of Lausanne. The project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie-Sklodowska Curie grant agreement No 838508. Maria Klepikova is grateful to Olivier Bour for valuable suggestions.

803

804

805

806

807
808 References

809 Anderson, M. P. (2005). Heat as a ground water tracer. *Ground water*, 810 43(6):951-968.

811 Andrade, J., Henrique, E., Almeida, M., and Costa, M. (2004). Heat trans-812 port through rough channels. *Physica A: Statistical Mechanics and its 813 Applications*, 339(3):296-310.

814 Bodvarsson, G. (1972). Thermal problems in the siting of reinjection wells. 815 *Geothermics*, 1(2):63 – 66.

816 Boschan, A., Auradou, H., Ippolito, I., Chertcoff, R., and Hulin, J. P. (2007). 817 Miscible displacement fronts of shear thinning fluids inside rough fractures. 818 *Water Resources Research*, 43(3).

819 Boschan, A., Ippolito, I., Chertcoff, R., Auradou, H., Talon, L., and Hulin, 820 J. P. (2008). Geometrical and taylor dispersion in a fracture with random 821 obstacles: An experimental study with fluids of different rheologies. *Water 822 Resources Research*, 44(6):W06420.

823 Bouchaud, E. (1997). Scaling properties of cracks. *Journal of Physics: Con-824 densed Matter*, 9(21):4319-4344.

825 Boutt, D. F., Grasselli, G., Fredrich, J. T., Cook, B. K., and Williams, 826 J. R. (2006). Trapping zones: The effect of fracture roughness on the 827 directional anisotropy of fluid flow and colloid transport in a single fracture. 828 *Geophysical Research Letters*, 33(21):L21402.

829 Brown, S. R. (1987). Fluid flow through rock joints: The effect of surface 830 roughness. *Journal of Geophysical Research: Solid Earth*, 92(B2):1337- 831 1347.

832 Brown, S. R., Kranz, R. L., and Bonner, B. P. (1986). Correlation be-833 tween the surfaces of natural rock joints. *Geophysical Research Letters*, 834 13(13):1430-1433.

835 Brush, D. J. and Thomson, N. R. (2003). Fluid flow in synthetic rough-walled 836 fractures: Navier-Stokes, Stokes, and local cubic law simulations. *Water 837 Resources Research*, 39(4).
838 Candela, T., Renard, F., Klinger, Y., Mair, K., Schmittbuhl, J., and Brodsky, E. E. (2012). Roughness of fault surfaces over nine decades of length scales. *Journal of Geophysical Research: Solid Earth*, 117(B8).

841 Cardenas, M. B., Slottke, D. T., Ketcham, R. A., and Sharp Jr., J. M. (2009). Effects of inertia and directionality on flow and transport in a rough asymmetric fracture. *Journal of Geophysical Research: Solid Earth*, 114(B6).

845 Chen, Y. and Zhao, Z. (2020). Heat transfer in a 3D rough rock fracture with heterogeneous apertures. *International Journal of Rock Mechanics and Mining Sciences*, 134:104445.

848 COMSOL Multiphysics User's Guide, Version 5.4. (2018). COMSOL AB, Stockholm, Sweden.

850 de Dreuzy, J.-R., Méheust, Y., and Pichot, G. (2012). Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (DFN). *J. Geophys. Res.: Solid Earth*, 117(B11):n/a–n/a.

854 de La Bernardie, J., de Dreuzy, J.-R., Bour, O., and Lesueur, H. (2019). Synthetic investigation of thermal storage capacities in crystalline bedrock through a regular fracture network as heat exchanger. *Geothermics*, 77:130–138.

858 Detwiler, R. L., Rajaram, H., and Glass, R. J. (2000). Solute transport in variable-aperture fractures: An investigation of the relative importance of Taylor dispersion and macrodispersion. *Water Resources Research*, 36(7):1611–1625.

862 Doetsch, J., Gischig, V., Krietsch, H., Villiger, L., Amann, F., Dutler, N., Jalali, R., Brixel, B., Klepikova, M., Roques, C., Giertzuch, P.-L., Kittilä, A., and Hochreutener, R. (2018). Grimsel ISC Experiment Description. Technical report, Zurich.

866 Fairley, J. P. (2009). Modeling fluid flow in a heterogeneous, fault-controlled hydrothermal system. *Geofluids*, 9(2):153–166.
868 Fox, D. B., Koch, D. L., and Tester, J. W. (2015). The effect of spatial aperture variations on the thermal performance of discretely fractured geothermal reservoirs. *Geothermal Energy*, 3(21).

871 Ge, S. (1997). A governing equation for fluid flow in rough fractures. *Water Resources Research*, 33(1):53–61.

873 Ge, S. (1998). Estimation of groundwater velocity in localized fracture zones from well temperature profiles. *Journal of Volcanology and Geothermal Research*, 84(1-2):93–101.

876 Geiger, S. and Emmanuel, S. (2010). W07504. *Water Resources Research*, 46(7):W07504.

878 Gossler, M. A., Bayer, P., and Zosseder, K. (2019). Experimental investigation of thermal retardation and local thermal non-equilibrium effects on heat transport in highly permeable, porous aquifers. *Journal of Hydrology*, 578:124097.

882 Gringarten, A. C., Witherspoon, P. A., and Ohnishi, Y. (1975). Theory of heat extraction from fractured hot dry rock. *Journal of Geophysical Research (1966-1977)*, 80(8):1120–1124.

885 Guo, B., Fu, P., Hao, Y., Peters, C. A., and Carrigan, C. R. (2016). Thermal drawdown-induced flow channeling in a single fracture in eggs. *Geothermics*, 61:62–62.

888 Hawkins, A. J., Fox, D. B., Becker, M. W., and Tester, J. W. (2017). Measurement and simulation of heat exchange in fractured bedrock using inert and thermally degrading tracers. *Water Resources Research*, 53(2):1210–1230.

892 He, X., Sinan, M., Kwak, H., and Hoteit, H. (2021). A corrected cubic law for single-phase laminar flow through rough-walled fractures. *Advances in Water Resources*, 154:103984.

895 Horne, R. N. and Rodriguez, F. (1983). Dispersion in tracer flow in fractured geothermal systems. *Geophysical Research Letters*, 10(4):289–292.

897 Incropera, F. P. and DeWitt, D. P. (1996). *Fundamentals of Heat and Mass Transfer*. John Wiley & Sons, Inc., New York City, New York, 4th edition edition.
900

901

902

903

904

Ishibashi, T., Watanabe, N., Hirano, N., Okamoto, A., and Tsuchiya, N. (2015). Beyond-laboratory-scale prediction for channeling flows through subsurface rock fractures with heterogeneous aperture distributions revealed by laboratory evaluation. *Journal of Geophysical Research: Solid Earth*, 120(1):106–124.

905

906

907

908

Johns, R. A., Steude, J. S., Castanier, L. M., and Roberts, P. V. (1993). Nondestructive measurements of fracture aperture in crystalline rock cores using x ray computed tomography. *Journal of Geophysical Research: Solid Earth*, 98(B2):1889–1900.

909

910

911

Jung, Y. and Pruess, K. (2012). A closed-form analytical solution for thermal single-well injection-withdrawal tests. *Water Resources Research*, 48(3):W03504.

912

913

914

915

Kant, M. A., Ammann, J., Rossi, E., Madonna, C., Höser, D., and Rudolf von Rohr, P. (2017). Thermal properties of central aare granite for temperatures up to 500 c: Irreversible changes due to thermal crack formation. *Geophysical Research Letters*, 44(2):771–776.

916

917

918

919

Klepikova, M. V., Borgne], T. L., Bour, O., Gallagher, K., Hochreutener, R., and Lavenant, N. (2014). Passive temperature tomography experiments to characterize transmissivity and connectivity of preferential flow paths in fractured media. *Journal of Hydrology*, 512:549 – 562.

920

921

922

923

Klepikova, M. V., Le Borgne, T., Bour, O., and Davy, P. (2011). A methodology for using borehole temperature-depth profiles under ambient, single and cross-borehole pumping conditions to estimate fracture hydraulic properties. *Journal of Hydrology*, 407(1):145–152.

924

925

926

927

928

Klepikova, M. V., Le Borgne, T., Bour, O., Dentz, M., Hochreutener, R., and Lavenant, N. (2016). Heat as a tracer for understanding transport processes in fractured media: Theory and field assessment from multiscale thermal push-pull tracer tests. *Water Resources Research*, 52(7):5442–5457.

929

930

Kocabas, I. (2005). Geothermal reservoir characterization via thermal injection backflow and interwell tracer testing. *Geothermics*, 34(1):27–46.
Kolditz, O. and Clauser, C. (1998). Numerical simulation of flow and heat transfer in fractured crystalline rocks: Application to the Hot Dry Rock site in Rosemanowes (U.K.). *Geothermics*, 27(1):1–23.

Lanahan, M. and Tabares-Velasco, P. C. (2017). Seasonal thermal-energy storage: A critical review on btes systems, modeling, and system design for higher system efficiency.

Lee, S. H., Yeo, I. W., Lee, K.-K., and Detwiler, R. L. (2015). Tail shortening with developing eddies in a rough-walled rock fracture. *Geophysical Research Letters*, 42(15):6340–6347.

Luo, J., Zhu, Y., Guo, Q., Tan, L., Zhuang, Y., Liu, M., Zhang, C., Xiang, W., and Rohn, J. (2017). Experimental investigation of the hydraulic and heat-transfer properties of artificially fractured granite. *Scientific Reports*, 7.

Luo, S., Zhao, Z., Peng, H., and Pu, H. (2016). The role of fracture surface roughness in macroscopic fluid flow and heat transfer in fractured rocks. *International Journal of Rock Mechanics and Mining Sciences*, 87:29–38.

Malkovsky, V. I. and Magri, F. (2016). Thermal convection of temperature-dependent viscous fluids within three-dimensional faulted geothermal systems: Estimation from linear and numerical analyses. *Water Resources Research*, 52(4):2855–2867.

Martinez, A. R., Roubinet, D., and Tartakovskiy, D. M. (2014). Analytical models of heat conduction in fractured rocks. *Journal of Geophysical Research: Solid Earth*, 119(1):83–98.

Méheust, Y. and Schmittbuhl, J. (2001). Geometrical heterogeneities and permeability anisotropy of rough fractures. *Journal of Geophysical Research: Solid Earth*, 106(B2):2089–2102.

Méheust, Y. and Schmittbuhl, J. (2003). Scale effects related to flow in rough fractures. *Pure and Applied Geophysics*, 160:1023–1050.

Méheust, Y. and Schmittbuhl, J. (2000). Flow enhancement of a rough fracture. *Geophysical Research Letters*, 27(18):2989–2992.
961

962

963

Neuville, A., Flekkøy, E. G., and Toussaint, R. (2013). Influence of asperities on fluid and thermal flow in a fracture: A coupled lattice boltzmann study. *Journal of Geophysical Research: Solid Earth*, 118(7):3394–3407.

964

965

966

967

968

Neuville, A., Toussaint, R., and Schmittbuhl, J. (2010a). Fracture roughness and thermal exchange: A case study at Soultz-sous-Forêts. *Comptes Rendus Geoscience*, 342(7):616 – 625. Vers l'exploitation des ressources géothermiques profondes des systèmes hydrothermaux convectifs en milieux naturellement fracturés.

969

970

Neuville, A., Toussaint, R., and Schmittbuhl, J. (2010b). Hydrothermal coupling in a self-affine rough fracture. *Phys. Rev. E*, 82:036317.

971

972

973

974

Neuville, A., Toussaint, R., and Schmittbuhl, J. (2011). Hydraulic transmissivity and heat exchange efficiency of open fractures: a model based on lowpass filtered apertures. *Geophysical Journal International*, 186(3):1064–1072.

975

976

977

Nick, H. M., Paluszny, A., Blunt, M. J., and Matthai, S. K. (2011). Role of geomechanically grown fractures on dispersive transport in heterogeneous geological formations. *Phys. Rev. E*, 84:056301.

978

979

Nottebohm, M., Licha, T., and Sauter, M. (2012). Tracer design for tracking thermal fronts in geothermal reservoirs. *Geothermics*, 43:37–44.

980

981

982

983

984

Odling, N. E. (1992). Permeability of natural and simulated fracture patterns. In Larsen, R., Brekke, H., Larsen, B., and Talleraas, E., editors, *Structural and Tectonic Modelling and its Application to Petroleum Geology*, volume 1 of *Norwegian Petroleum Society Special Publications*, pages 365–380. Elsevier, Amsterdam.

985

986

Oron, A. P. and Berkowitz, B. (1998). Flow in rock fractures: The local cubic law assumption reexamined. *Water Resources Research*, 34(11):2811–2825.

987

988

989

990

Pandey, S., Chaudhuri, A., Kelkar, S., Sandeep, V., and Rajaram, H. (2014). Investigation of permeability alteration of fractured limestone reservoir due to geothermal heat extraction using three-dimensional thermo-hydrochemical (THC) model. *Geothermics*, 51(Supplement C):46 – 62.
991

992

993

Patterson, J. W. and Driesner, T. (2020). The effect of thermo-elastic stress re-distribution on geothermal production from a vertical fracture zone. *Geothermics*, 85:101745.

994

995

996

Patterson, J. W., Driesner, T., Matthai, S., and Tomlinson, R. (2018). Heat and fluid transport induced by convective fluid circulation within a fracture or fault. *Journal of Geophysical Research: Solid Earth*, 123(4):2658–2673.

997

998

999

Plouraboué, F., Kurowski, P., Boffa, J.-M., Hulin, J.-P., and Roux, S. (2000). Experimental study of the transport properties of rough self-affine fractures. *Journal of Contaminant Hydrology*, 46(3):295–318.

1000

1001

1002

Pruess, K. and Doughty, C. (2010). Thermal single-well injection-withdrawal tracer tests for determining fracture-matrix heat transfer area. *Conference: 35th Workshop on Geothermal Reservoir Engineering, Stanford, CA, USA*.

1003

1004

1005

1006

1007

Read, T., Bour, O., Bense, V., Le Borgne, T., Goderniaux, P., Klepikova, M. V., Hochreutener, R., Lavenant, N., and Boschero, V. (2013). Characterizing groundwater flow and heat transport in fractured rock using fiber-optic distributed temperature sensing. *Geophysical Research Letters*, 40(10):2055–2059.

1008

1009

1010

Renard, F., Candela, T., and Bouchaud, E. (2013). Constant dimensionality of fault roughness from the scale of micro-fractures to the scale of continents. *Geophysical Research Letters*, 40(1):83–87.

1011

1012

1013

Saar, M. O. (2010). Review: Geothermal heat as a tracer of large-scale groundwater flow and as a means to determine permeability fields. *Hydrogeology Journal*, 19(1):31–52.

1014

1015

1016

Salimzadeh, S. and Nick, H. (2019). A coupled model for reactive flow through deformable fractures in enhanced geothermal systems. *Geothermics*, 81:88–100.

1017

1018

1019

Salimzadeh, S., Paluszny, A., Nick, H. M., and Zimmerman, R. W. (2018). A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems. *Geothermics*, 71:212–224.

1020

1021

Schmittbuhl, J., Gentier, S., and Roux, S. (1993). Field measurements of the roughness of fault surfaces. *Geophysical Research Letters*, 20(8):639–641.
Shortall, R., Davidsdottir, B., and Axelsson, G. (2015). Geothermal energy for sustainable development: A review of sustainability impacts and assessment frameworks. *Renewable and Sustainable Energy Reviews*, 44(C):391–406.

Silliman, S. and Robinson, R. (1989). Identifying fracture interconnections between boreholes using natural temperature profiling: I. conceptual basis. *Groundwater*, 27(3):393–402.

Taron, J. and Elsworth, D. (2009). Thermal-hydrologic-mechanical-chemical processes in the evolution of engineered geothermal reservoirs. *International Journal of Rock Mechanics and Mining Sciences*, 46(5):855–864.

Tsang, C.-F. (1991). Coupled hydromechanical-thermochemical processes in rock fractures. *Reviews of Geophysics*, 29(4):537–551.

Tsang, Y. W. (1992). Usage of “equivalent apertures” for rock fractures as derived from hydraulic and tracer tests. *Water Resources Research*, 28(5):1451–1455.

Tsang, Y. W. and Tsang, C. F. (1987). Channel model of flow through fractured media. *Water Resources Research*, 23(3):467–479.

Vik, H., Salimzadeh, S., and Nick, H. M. (2018). Heat recovery from multiple-fracture enhanced geothermal systems: The effect of thermoelastic fracture interactions. *Renewable Energy*, 121:606–622.

Wang, L. and Cardenas, M. B. (2014). Non-fickian transport through two-dimensional rough fractures: Assessment and prediction. *Water Resources Research*, 50(2):871–884.

Wang, L. and Cardenas, M. B. (2017). Transition from non-fickian to fickian longitudinal transport through 3-d rough fractures: Scale-(in)sensitivity and roughness dependence. *Journal of Contaminant Hydrology*, 198:1–10.

Wang, L., Cardenas, M. B., Slottke, D. T., Ketcham, R. A., and Sharp Jr., J. M. (2015). Modification of the local cubic law of fracture flow for weak inertia, tortuosity, and roughness. *Water Resources Research*, 51(4):2064–2080.
1054

1055

1056

1057

Wang, L., Cardenas, M. B., Zhou, J.-Q., and Ketcham, R. A. (2020). The complexity of nonlinear flow and non-fickian transport in fractures driven by three-dimensional recirculation zones. *Journal of Geophysical Research: Solid Earth*, 125(9):e2020JB020028. e2020JB020028 2020JB020028.

1058

1059

1060

Yang, Z., Méheust, Y., Neuweiler, I., Hu, R., Niemi, A., and Chen, Y.-F. (2019). Modeling immiscible two-phase flow in rough fractures from capillary to viscous fingering. *Water Resources Research*, 55(3):2033–2056.

1061

1062

Yoon, S. and Kang, P. K. (2021). Roughness, inertia, and diffusion effects on anomalous transport in rough channel flows. *Phys. Rev. Fluids*, 6:014502.

1063

1064

1065

Zhang, C.-L., Conil, N., and Armand, G. (2017). Thermal effects on clay rocks for deep disposal of high-level radioactive waste. *Journal of Rock Mechanics and Geotechnical Engineering*, 9(3):463–478.

1066

1067

1068

1069

Zhou, Q., Oldenburg, C. M., Rutqvist, J., and Birkholzer, J. T. (2017). Revisiting the fundamental analytical solutions of heat and mass transfer: The kernel of multirate and multidimensional diffusion. *Water Resources Research*, 53(11):9960–9979.

1070

1071

1072

Zimmerman, R. W. and Yeo, I.-W. (2013). *Fluid Flow in Rock Fractures: From the Navier-Stokes Equations to the Cubic Law*, pages 213–224. American Geophysical Union (AGU).
Author statement

**MK:** Conceptualization, Methodology, Software, Writing- Original draft preparation, Writing- Reviewing and Editing, Funding acquisition

**YM:** Conceptualization, Methodology, Software, Editing

**NL:** Conceptualization, Editing, Supervision

**CR:** Methodology, Software

<figure><img src="image_15.png" /></figure>
**Declaration of interests**

[x] The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

[ ] The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

<figure><img src="image_16.png" /></figure>