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

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 lubri- fication 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. fracturing process itself, as well as post-fracturing processes such as geological stress and strain, chemical dissolution, precipitation and


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; 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 (de La Bernardie et al., 2019;Lanahan and Tabares-Velasco, 2017) and near-field thermal effects in the context of radioactive waste disposal . 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 and Tsang, 1987;Méheust and Schmittbuhl, 2000;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 Boutt et al., 2006;Cardenas et al., 2009;Ge, 1997;Méheust and Schmittbuhl, 2001, Wang and Cardenas, 2014Yang 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 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 m is shown by the dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

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). Parameter Matrix, r Water, f Density ρ, kg/m 3 2500 1000 Thermal conductivity k, W/(m K) 0.59 3.5 Heat capacity Cp, J/(kg K) 750 4189 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. (2010bNeuville et al. ( , 2010aNeuville et al. ( , 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 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, γ/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, γ/d m = 0.6 (c); for the fracture demonstrating the smallest ratio d h /d m = 0.18 at the largest closure considered, γ/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). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 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 and in a model with strong fracture roughness belonging to fracture family A (see Fig. 2c), with γ/d m = 0.46 and L/Lc = 1 ((b) -fracture plane view and (d) -mid-longitudinal crosssection 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 flush = 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). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 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 ζ = 0.8d m , with aperture closures γ/d m varying from 0.001 to 0.6 over a wide range of mean fracture apertures d m , 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;2010a;, 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.

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, 2018), as detailed in the following.

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 ≤ 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 γ of the aperture field over the entire fracture. In our model, the Hurst (or, roughness) exponent has been chosen constant at ζ = 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 γ/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 Fig. 1b.

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 . 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, 2000): where P is the local pressure [Pa] and η 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 (ρ f l 2 z )/(ηl h ) (Méheust and Schmittbuhl, 2001;Neuville et al., 2011), where ρ f is the fluid density [kg/m 3 ], 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 = ∇Pd 2 m /8η 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 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 (∇⋅u = 0), which implies that q is also conservative: ∇⋅q = 0 (Plouraboué et al., 2000). Inserting the local cubic law (Eq. 1) in this conservation law yields the Reynolds equation: 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 ∇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. Odling, 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 Fig. 1b, is shown in Fig. 1a as black arrows.
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 fracture 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 Δ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.

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. Fig. 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, 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: where the conductive heat flux is given by In the fracture: where the conductive heat flux is given by Here ∇ t denotes the gradient operator restricted to the fracture's tangential plane, ρ the density [kg/m 3 ], C p the heat capacity [J/kgK], k the thermal conductivity [W/mK], 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): 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 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 lubri-cation 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;. Here and below, thermal parameters are selected to represent granite, the host formations of most deep geothermal projects (Table 1).

Hydraulic behaviour
In Fig. 2a we present the ratio of hydraulic to mechanical apertures d h /d m , as a function of the fracture closure γ/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 γ/d m and/or values for the mismatch scale L c (investigated in Section 3.2.5). In Fig. 2a Méheust and Schmittbuhl, 2001;, 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 Figs. 2b-d some of the investigated aperture fields d fr (x, y)/d m for the highest closure considered, and in Figs. 2eg the corresponding maps of local fluxes (2-D velocities) normalized by their mean value (q/ ‖ q ‖). In Fig. 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 Fig. 2e. This fracture morphology corresponds to the largest ratio d h /d m in Fig. 2a. Fig. 2c shows the map of the ratio d fr (x, y)/d m for one of the fracture families considered in Fig. 2a; this family demonstrates a moderate flowinhibiting behavior (family A, red markers). In this case, the fracture aperture field is characterized by a main tortuous channel with smaller flow obstacles (Fig. 2f). In Fig. 2d a barrier is seen across the whole fracture that is perpendicular to the applied pressure gradient. As shown in Fig. 2g, this results in a strong flow-inhibiting behaviour of the fracture (lowest ratio d h /d m in Fig. 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). Fig. 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 Fig. 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

Thermal front definition
The temperature distribution in the rough-walled fracture is highly heterogeneous (Fig. 3b) and the temperature evolution over time may differ considerably even for points located close to each other (Fig. 3e). Initially, the thermal anomaly propagates along preferential large aperture channels and reaches for instance points B, D and C of Fig. 3b, whose temperature evolution is shown in Fig. 3e. 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 Fig. 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 Fig. 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 Fig. 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 Fig. 2 for various closure γ/d m values, ζ = 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. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 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 ΔT = 1 /2 (black lines in Fig. 3).
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 Fig. 3) that is delayed relative to the fluid front (blue solid line in Fig. 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 (Fig. 3a).
In order to characterize how the thermal behavior evolves on average, we consider the evolution of the thermal front with time. Fig. 3b presents an example of the thermal front advancement for the fracture family A, γ/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.

Influence of fracture closure
We now investigate how the roughness amplitude influences the thermal behaviour of fractures from family A (Fig. 2c), which exhibits a moderate flow-inhibiting behavior (red markers in Fig. 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 Fig. 4a. For small fracture closure γ/d m = 0.02 (Fig. 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 (Fig. 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 Fig. 4b, where maps of local fluxes normalized by the mean local flux are shown. In Fig. 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 Fig. 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 Fig. 2. Finally, Fig. 4c presents the simulated temperature fields (after 1000 s of injection) for different values of fracture closure.
As discussed above, one of the important characteristics of the ge- ometry of the mixing zone (where 0 < ΔT < 1) is the position and shape of the thermal front (black dashed lines in Fig. 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 γ/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 γ/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 (Fig. 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)   front position, which quantifies its width along the mean flow direction, is plotted in Fig. 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 Fig. 2a.
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). Fig. 6 presents the evolution in time of the average velocity of the thermal front v ROUGH = x * ROUGH /t relative to the v PP = x * PP /t, that is, the thermal front velocity in a fracture modeled with a constant aperture d m . For a small fracture closure γ/d m = 0.05, which corresponds to a nearly smooth aperture field, a parallel plate model reproduces a similar thermal profile and the ratio v ROUGH /v PP is close to 1. As the fracture closure is increased, γ/d m = 0.18, 0.32, 0.46 and 0.60 (Fig. 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 .
For all simulations, higher ratios of thermal front velocities v ROUGH / v PP are observed close to the fracture inlet. Fig. 6 also demonstrates that for small fracture closure γ/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 γ /d m (Fig. 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, γ/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 Fig. 6) with large apertures (d m = 13 − 23 mm), the thermal front velocity ratio becomes quasi-steady after t = 400 − 700 s, respectively.

Influence of the fracture hydraulic aperture
For the same cases as in Fig. 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 in the channeled areas of a rough-walled fracture.
However, once the hosting rock temperature starts to evolve, a process which was not accounted for in the model of Neuville et al. (2010b), the thermal front velocity in rough fractures slows down (Fig. 6), and, for t > 400 s, we obtain (Fig. 7) a perfect correlation between the ratio v ROUGH /v PP and the ratio between the hydraulic and mechanical apertures d h /d m (black triangle markers). This implies that once heat has diffused along the out-of-plane direction over the fracture's aperture, the thermal behavior of a rough-walled fracture is determined by its effective hydraulic transmissivity. For the morphology investigated here (Family A fractures), the permeability is reduced, and, thus, the advancement of the thermal front is slower in rough fractures compared to what would be expected with a model of flat fractures having the same mechanical aperture d m (v ROUGH /v PP < 1).
We now leave the specifics of family A and consider all fracture families in Fig. 2a for different closure values (Test 7, Table 2). Statistical thermal results presented in Fig. 8 confirm the near-perfect correlation between the ratio v ROUGH /v PP and the ratio between the hydraulic and mechanical apertures d h /d m . 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 roughwalled fracture is determined by the hydraulic aperture d h . For different fracture families with the same fracture closure γ/d m , we find that the flow inhibiting behavior 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 ROUGH /v PP < 1.6) from one fracture family to another, but ratios smaller than 1 are favored statistically.

Conductive heat flux
While previous studies characterized the heat exchange efficiency of rough fractures through the use of temperature metrics (Neuville et al., 2010a;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, γ/d m = 0.05, 0.18, 0.32, 0.46, 0.60. Our results, presented in Fig. 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 PP converges in time to a plateau (Fig. 9, inset). Interestingly, Fig. 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 Table 3 Péclet numbers considered in previous studies on heat transport in fractured media.

Study
Scale Peclet number Horne and Rodriguez (1983) Fracture network 10-100 Ge (1998) Single fracture 5-100 Geiger and Emmanuel (2010) Fracture network 14-145 Neuville et al. (2010a) Single fracture 7 000 Neuville et al. (2013) Single fracture 10 Klepikova et al. (2016) Single fracture 5 000 Hawkins et al. (2017) Single fracture 2-8 de La Bernardie et al. (2019) Single fracture 70 000 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 (Fig. 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).

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 γ/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). Fig. 10 presents the simulated temperature fields (after 1000 s of injection) for different mismatch scales, L /L c = 2 (Fig. 10a), L/L c = 4 (Fig. 10b), L/L c = 16 (Fig. 10c), and for different values of fracture closure, γ/d m = 0. 02, 0.25, 0.48, and 0.60 (from top to bottom). Similarly to L/L c = 1 (Fig. 4), the thermal channeling effect is all the more pronounced as the fracture is more closed (Fig. 10). For L / L c = 2 (Fig. 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 (Fig. 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.
Despite the changes in the q-field with increasing L /L c , we find similarly to what is observed for L c = L (black markers in Fig. 8), that the fluid flow scaling translates into that of heat transport. Indeed, red markers in Fig. 8 show strong correlation between the ratio v ROUGH /v PP and the ratio between the hydraulic and mechanical apertures d h /d m for L/L c = 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 c .

Discussion
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).
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;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 temperaturedependence 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 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 acounting 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).
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.
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.

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 ζ = 0.8, with aperture closures γ /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 displacement of the thermal front along the fracture and less thermal exchange between the fracture and the surrounding rock.
Our main finding is that under the considered conditions, thermal behaviour of rough-walled rock fractures only depends on the hydraulic properties. 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 fractures 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 transition 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 ∼ 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 aperture and boundary conditions. This striking novel finding results from an improved description of the coupled flow and heat transport.
The practical implication of our finding is that thermal exchanges at the scale of a single fracture is controlled by the effective hydraulic transmissivity. Provided that thermal properties of the host rock are known, this implies that (1) heat tracer tests are reliable for inferring effective fracture transmissivity, 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 predicted over significantly large time scales using well-known low-dimensional hydraulic parameterizations in terms of effective hydraulic properties.
Future work could address non-Stokes flow conditions (i.e., laminar but non-linear flow) in the fracture. Another interesting prospect is the overall 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. The study of coupled flow and heat transport in such DFNs will be the topic of future work.