## 1 Introduction

Modeling of the fluid transport in low permeable crustal rocks is of central importance for many applications (Neuman, 2005). Among them is the monitoring of the geothermal circulation in the project of Soultz-sous-Forêts, France, (Bachler et al., 2003) where the heat exchange especially occurs through open fractures in granite (Gérard et al., 2006).

Numerous hydrothermal models have already been proposed. For simple geometries, some analytical solutions are known: e.g., the cases of parallel plates (Turcotte and Schubert, 2002) or flat cylinders (Heuer et al., 1991). More complex models exist as well like the models of three-dimensional (3D) networks of fractures reproducing geological observations and possibly completed with stochastical distributions of fractures (e.g. in Soultz-sous-Forêts, France, (Gentier, 2005; Rachez et al., 2007) or in Rosemanowes, UK (Kolditz and Clauser, 1998)). Nevertheless, the geometry of each fracture is generally simple. Kolditz and Clauser (1998) have however suspected that differences between heat models and field observations could be due to channeling induced by the fracture roughness or the fracture network. Channeling of the fluid flow owing to fracture roughness has indeed already been experimentally observed and studied (Méheust and Schmittbuhl, 2000; Plouraboué et al., 2000; Schmittbuhl et al., 2008; Tsang and Tsang, 1998).

Here, we limit our study to the fracture scale and we will show only one example of thermal behavior, among other simulations we completed (Neuville et al., submitted). The specificity of our hydrothermal model is to take into account the different scale fluctuations of the fracture morphology. We aim at bringing out the main parameters, which control the hydraulic and thermal behavior of a complex rough fracture. The perspective is to propose a small set of effective parameters that could be introduced within simplified elements for an upscaled network model.

We first describe our geometrical model of the fracture aperture thanks to self-affine apertures. Then, using lubrication approximations, we obtain the bidimensional (2D) pressure and thermal equations when a cold fluid is injected through the fracture in a stationary regime. The temperature within the surrounding rock is supposed to be hot and constant in time and space. The fluid density is also supposed to be constant.

We apply our numerical model to the case study at Soultz-sous-Forêts and we show for this case an example of the computed hydraulic and thermal behavior. Finally, we aim at bringing out what is the minimal geometrical information needed to get the dominant behavior of the hydraulic and thermal fields. This last approach is based on spatial low pass Fourier filtering of the geometrical aperture field.

## 2 Modeling

### 2.1 Roughness of the fracture aperture

We consider that the mean fracture plane is described by the $(\stackrel{\u02c6}{x},\stackrel{\u02c6}{z})$ coordinates and the perpendicular direction is $\stackrel{\u02c6}{y}$ (Fig. 1) – where the hat notation refers to unit vectors along the (x,y,z) axis. It has been shown that a possible geometrical model of natural rough fractures consists in self-affine surfaces. A surface described by a function y = f(x,z) is self-affine if it is statistically invariant under the scaling transformation x → λx, z → λz and y → λ^{ζ}y, where ζ is called the roughness exponent or Hurst exponent. Such surfaces are therefore statistically invariant upon an isotropic scaling within their mean plane while along the perpendicular direction, the scaling is anisotropic (e.g. Brown and Scholz, 1985; Cox and Wang, 1993; Power et al., 1987; Schmittbuhl et al., 1993, 1995). Most fracture surfaces in heterogeneous material exhibit a Hurst exponent equal to ζ = 0.8 (Bouchaud, 1997; Santucci et al., 2007; Schmittbuhl et al., 1993, 1995). Sandstone fractures, however, show ζ = 0.5 (Boffa et al., 1998; Méheust, 2002).

It is important to note that a self-affine surface having a roughness exponent smaller than one is asymptotically flat at large scales (Roux et al., 1993). Accordingly, the self-affine topography can be seen as a perturbation of a flat interface. When the lubrication approximation (Pinkus and Sternlicht, 1961) holds, in particular with smooth enough self-affine perturbations or highly viscous fluid, only the local aperture controls the flow and not the local slope of the fracture. The accuracy of the lubrication approximation, compared to the full Navier-Stokes resolution, was studied in Al-Yaarubi et al. (2005). Under this assumption, the only required geometrical input is the aperture field (also called the geometrical aperture); there is no need to know the geometry of each facing fracture surfaces. The aperture between two uncorrelated self-affine fracture surfaces having the same roughness exponent is as well self-affine (Méheust and Schmittbuhl, 2003). Thus, we generate the numerical apertures by using self-affine functions.

Several independent self-affine aperture morphologies can be generated with the same roughness exponent chosen equal to ζ = 0.8. They exhibit various morphology patterns according to the chosen seed of the random generator (Méheust, 2002). The mean geometrical aperture A and the root-mean square deviation σ (RMS) of an aperture a(x,z) are defined as

$$A=\frac{\iint a\text{\hspace{0.28em}}dxdz}{{l}_{x}\cdot {l}_{z}}\text{and}$$ | (1) |

$${\sigma}^{2}=\frac{\iint {a}^{2}dxdz}{{l}_{x}\cdot {l}_{z}}-{A}^{2}$$ | (2) |

_{x}the length and l

_{z}the width of the fracture. To keep the boundary geometry of the domain as simple as possible, we do not allow any contact area (i.e. no local aperture equal to zero). This is obtained by considering a large enough aperture average to get strictly positive aperture fields.

It has to be noted that our hydrothermal model can be applied to other geometrical models (i.e. different from a self-affine model), which might be more relevant depending on the geological context.

### 2.2 Physics of hydraulic flow

The hydraulic flow is obtained under the same hypotheses and solved in the same way as in Méheust and Schmittbuhl (2001). We use finite differences, and the system of linear equations is inverted using an iterative biconjugate gradient method (Press et al., 1992).

We impose a pressure drop across the system and study the steady state flow of a Newtonian fluid at low Reynolds number, so that the viscous term dominates the inertial one in the Navier-Stokes equation (Batchelor, 2002; Stokes, 1846):

$$\overrightarrow{\nabla}P=\eta \Delta {\overrightarrow{u}}^{3D}$$ | (3) |

^{3D}the velocity of the fluid and P is the pressure deviation from the hydrostatic profile (or the hydraulic head equal to the pressure corrected by the gravity effect). To be in the framework of the lubrication approximation (Pinkus and Sternlicht, 1961), we consider fractures with constant enough apertures together with a small Reynolds number. In doing so, the velocity vector of the fluid flow has negligible components normal to the mean fracture plane. We consider that the macroscopic pressure gradient is imposed along $\stackrel{\u02c6}{x}$; $\stackrel{\u02c6}{z}$ is therefore perpendicular to the mean flow direction. Accordingly, the fluid velocity follows a parabolic law (e.g. Iwai, 1976) (Fig. 2):

$${\overrightarrow{u}}^{3D}(x,y,z)=\frac{{\overrightarrow{\nabla}}_{2}P}{2\eta}(y-{y}_{1})(y-{y}_{2})$$ | (4) |

_{1}and y

_{2}are the local fracture sides coordinates and ${\overrightarrow{\nabla}}_{2}$ is the gradient operator in the fracture plane. The hydraulic flow through the fracture aperture follows a cubic law:

$$\overrightarrow{q}(x,z)={\int}_{a}{\overrightarrow{u}}^{3D}(x,y,z)dy=-\frac{{a}^{3}}{12\eta}{\overrightarrow{\nabla}}_{2}P$$ | (5) |

and the bidimensional (2D) velocity $\overrightarrow{u}$ is defined from the average of the velocity ${\overrightarrow{u}}^{3D}$ over the aperture with

$$\overrightarrow{u}(x,z)=\frac{1}{a(x,z)}{\int}_{a}{\overrightarrow{u}}^{3D}(x,y,z)dy=-\frac{{a}^{2}}{12\eta}{\overrightarrow{\nabla}}_{2}P$$ | (6) |

Furthermore, considering the fluid to be incompressible, the Reynolds equation is obtained: ${\overrightarrow{\nabla}}_{2}({a}^{3}{\overrightarrow{\nabla}}_{2}P)=0$. As boundary conditions of this equation, we impose the pressure at the inlet and outlet of the fracture (if x = 0, P = P_{0} and if x = l_{x}, P = P_{lx}, with P_{0} > P_{lx}) and consider impermeable sides at z = 0 and z = l_{z}.

### 2.3 Physics of thermal exchange

On the basis of a classical description (e.g. Ge, 1998; Turcotte and Schubert, 2002), we aim at modeling the fluid temperature when cold water is permanently injected at the inlet of a hot fracture at temperature T_{0}. As the conduction inside the rock is not taken into account (hypothesis of infinite thermal conduction inside the rock), the fracture sides are supposed to be permanently hot at the fixed temperature T_{r}. This hypothesis should hold for moderate time scales (e.g. minutes), after the fluid injection has stabilized, and before the rock temperature has significantly changed, or alternatively once the whole temperature bedrock is stabilized (which depends on the boundary condition of the entire region). For time scales implying evolution of the rock temperature, our model should be coupled to a model of the rock temperature evolution.

The fluid temperature is controlled by the balance between thermal convection and conduction inside the fluid, which reads (Landau and Lifchitz, 1994): ${\overrightarrow{u}}^{3D}\cdot \overrightarrow{\nabla}T=\chi \Delta T$ where χ is the thermal diffusivity of the fluid and T the fluid temperature. We extend the local lubrication approximation by considering that the slopes of the fracture morphology are small enough to limit the conduction only along the y-axis. We suppose that the leading terms are the conduction along the y-axis and the in-plane convection (since there is no fluid velocity component along $\stackrel{\u02c6}{y}$). Indeed, the off-plane free convection has been shown to be negligible (its magnitude is of the order of km/year (Bataillé et al., 2006)). So, the previous equation reduces to:

$$\frac{{\partial}^{2}T}{\partial {y}^{2}}=\frac{{\overrightarrow{u}}_{x}^{3D}}{\chi}\frac{\partial T}{\partial x}+\frac{{\overrightarrow{u}}_{z}^{3D}}{\chi}\frac{\partial T}{\partial z}$$ | (7) |

_{x}(∂T/∂x)+q

_{z}(∂T/∂x) is independent of y, where q

_{x}and q

_{z}are the in-plane component of $\overrightarrow{q}$ defined in Eq. (5). Accordingly, we find that the temperature solution has a quartic profile (Fig. 2) along the fracture aperture:

^{1}

$$T=-\frac{\beta}{2{a}^{3}\chi}(y-{y}_{1})(y-{y}_{2})(y-\sqrt{5}{y}_{1})(y-\sqrt{5}{y}_{2})+{T}_{r}$$ | (8) |

_{1}and y

_{2}are the local fracture sides coordinates.

Similarly to what is done for the hydraulic flow, we solve the thermal equation by integrating it along the fracture aperture (following the lubrication approximation extended to the thermal field). In particular, when doing the balance of the energy fluxes, we express the advected free energy flux as $\rho c{\int}_{a}{u}^{3D}(x,y,z)[T(x,y,z)-{T}_{0}]dy$. Accordingly, we introduce:

$$\overline{T}(x,z)=\frac{{\int}_{a}{u}^{3D}(x,y,z)T(x,y,z)dy}{{\int}_{a}{u}^{3D}(x,y,z)dy}$$ | (9) |

_{r}/φ

_{ref}which compares the efficiency of the heat flow along the fracture boundaries:

$${\phi}_{r}=\mp \chi \rho c{\left.\frac{\partial T}{\partial y}\right|}_{y={y}_{1},{y}_{2}}$$ | (10) |

$$\overrightarrow{q}\cdot {\overrightarrow{\nabla}}_{2}\overline{T}+2\frac{\chi}{a}Nu(\overline{T}-{T}_{r})=0\text{,}$$ | (11) |

Boundary conditions are: $\overline{T}(0,z)={T}_{0}$ at the inlet and $\overline{T}({l}_{x},z)={T}_{r}$ at the outlet (with l_{x} large enough). Any boundary condition for the temperature along z = 0 or z = l_{z} can be used as the hydraulic flow $\overrightarrow{q}$ is null there.

We discretize this equation by using a first order finite difference scheme and finally get $\overline{T}$ by inverting the system using a biconjugated gradient method (Press et al., 1992).

It is finally possible to get the three-dimensional temperature field T anywhere within the fluid by using the previous β expression and the quartic profile (Eq. (8). Fig. 3 illustrates an example of temperature field at a given z = z_{0}, T(x,y,z = z_{0}), obtained in that way from a given bidimensional field $\overline{T}$. Along any given cut at x = x_{0}, the temperature (represented by the color scale) follows a quartic law (Fig. 2). The boundaries between the colors are isotherms.

### 2.4 Definition of characteristic quantities describing the computed hydraulic and thermal fields

#### 2.4.1 Comparison to modeling without roughness

If we consider a fracture modeled by two parallel plates separated by a constant aperture A, then the gradient of pressure is constant all along the fracture as well as the hydraulic flow which is equal to:

$${\overrightarrow{q}}_{//}=-\frac{\Delta P}{{l}_{x}}\frac{{A}^{3}}{12\eta}\stackrel{\u02c6}{x}$$ | (12) |

_{lx}− P

_{0}. Under these conditions, the analytical solution of Eq. (11) is:

$${\overline{T}}_{//}=({T}_{0}-{T}_{r})\mathrm{exp}\left(-\frac{x}{{R}_{//}}\right)+{T}_{r}$$ | (13) |

_{//}is a thermal length describing the distance at which the fluid typically reaches the temperature of the surrounding rock. We have:

$${R}_{//}=\frac{A\left|\right|{\overrightarrow{q}}_{//}\left|\right|}{2Nu\chi}=-\frac{\Delta P}{{l}_{x}}\frac{{A}^{4}}{24\eta Nu\chi}=\frac{APe}{Nu}$$ | (14) |

For rough fractures, we want to study whether the temperature profiles along x at a coarse grained scale can still be described by Eq. (13) and, if so, what is the impact of the fracture roughness on the thermal length R.

#### 2.4.2 Hydraulic aperture

The hydraulic flow can be macroscopically described using the hydraulic aperture H (Brown, 1987; Zimmerman et al., 1991), defined as the equivalent parallel plate aperture to get the macroscopic flow 〈q_{x}〉 under the pressure gradient ΔP/l_{x}:

$$H={\u2329-{q}_{x}\frac{12\eta {l}_{x}}{\Delta P}\u232a}^{1/3}\text{,}$$ | (15) |

$$h(x,z)={\left(-{q}_{x}(x,z)\frac{12\eta {l}_{x}}{\Delta P}\right)}^{1/3}$$ | (16) |

Local geometrical and hydraulic apertures are denoted here with small letters while the corresponding macroscopic variables (mean geometrical and hydraulic aperture) are in capital letters.

#### 2.4.3 Thermal aperture

For the thermal aspect, once $\overline{T}$ is known, we aim at defining a thermal length R like in Eq. (13). To do that, we define $\overline{\overline{T}}$, a z-average temperature which varies only along the forced gradient direction x, and weighted by the 2D fluid velocity u_{x} to fulfill energy conservation:

$$\overline{\overline{T}}(x)=\frac{\underset{{l}_{z}}{\int}{u}_{x}(x,z)\overline{T}(x,z)dz}{{\int}_{{l}_{z}}{u}_{x}(x,z)dz}$$ | (17) |

Then, based on the flat plate temperature solution (Eq. (13), we do a linear fit of $ln[(\overline{\overline{T}}-{T}_{r})/({T}_{0}-{T}_{r})]$ plotted as a function of x, and we use the slope of this fit to get the characteristic thermal length R. This fit is computed over the zone where the numerical precision of the fitted quantities is sufficient (larger than ln[2 × 10^{−6}]). Similarly to the parallel plate case (Eq. (14)), the thermal length R can be used to define a thermal aperture Γ:

$$R=-\frac{\Delta P}{{l}_{x}}\frac{{\Gamma}^{4}}{24\eta Nu\chi}$$ | (18) |

which means that a fracture modeled by parallel plates separated by a distance Γ provides the same averaged thermal behavior as the rough fracture of mean geometrical aperture A.

## 3 Case study at Soultz-sous-Forêts (France)

### 3.1 Computation of apertures, hydraulic and thermal fields

Let us consider the GPK3 and GPK2 wells of the deep geothermal drilling near Soultz-sous-Forêts (France), which are separated by a distance of about 600 m at roughly 5000 m of depth. From hydraulic tests (Sanjuan et al., 2006), it has been shown that the hydraulic connection between both wells is relatively direct and straight. Sausse et al. (Sausse et al., 2008) showed that actually a fault is linking GPK3 (at 4775 m) to GPK2. This fault zone consists of a large number of clusters of small fractures (most apertures ranges around 0.1 mm) (Sausse et al., 2008) which probably lead to complex hydraulic streamlines and heat exchanges. We study here a simplified model of this connecting fault zone between the wells using one single rectangular rough fracture. The size of the studied fracture is l_{x} × l_{y} = 680 × 1370 m^{2}. Individual fracture apertures are typically of the order of 0.2 mm (Genter and Jung, private communication) while the fracture zone is rather thicker (10 cm) (Sausse et al., 2008). To account for this variability of the fault zone aperture, we use a probabilistic model with the following macroscopic properties: a mean aperture A equal to 3.60 mm and its standard deviation to σ = 1.23 mm. Fig. 4 shows an example of a self-affine aperture randomly generated with the required parameters.

With little knowledge about the pressure conditions along the boundaries of this model, we assume that the two facing sides along x of this rectangular fracture correspond to the inlet and outlet of the model where the pressure is homogeneous, respectively P_{0} and P_{lx}. In other words, we assume the streamlines to be as straight as possible between both wells.

The pressure gradient ΔP/l_{x} is chosen as −10^{−2} bar/m, which corresponds to about six bars between the bottom of both wells. The dynamic viscosity is chosen to be 3 × 10^{−4} Pa/s (reference value for pure water at 10 Pa and 100 °C in Spurk and Aksel (2008)). The Reynolds number rescaled with the characteristic dimensions of the fracture (Méheust and Schmittbuhl, 2001) is equal to Re′ = (ρu_{x} a^{2})/(η.l_{x}) = 0.026 and the Péclet number is Pe = 3.8 × 10^{4}.

Then, we solve the hydraulic flow in the fracture domain and obtain the 2D velocity field, $\overrightarrow{u}$ defined in Eq. (6). Fig. 5 shows the spatial fluctuations of $\left|\right|\overrightarrow{u}\left|\right|$. For information, a parallel plate model separated by the chosen aperture A would predict a homogeneous fluid velocity of 3.6 m/s and a thermal length R_{//} = 33.3 m. As we see in Fig. 5, the 2D velocity field exhibits interesting features: the fluid is rather immobile along the upper and lower borders of the fracture (close to z = 0 and z = l_{x}) while most of the fluid flows very quickly through a channel in the middle of the fracture.

The macroscopic hydraulic aperture is deduced from the local hydraulic flow estimate (Eq. (15)): H = 3.73 mm, which is slightly higher than the mean mechanical aperture A = 3.60 mm. Therefore, this fracture is more permeable than parallel plates separated by A. In other words, the fracture is geometrically thinner than what one would expect from the knowledge of H possibly inverted from an hydraulic test. However, the local hydraulic apertures h (Eq. (16)) range from nearly 0 to 5.43 mm (Fig. 6).

From the average estimate of the fluid velocity, we can go back to our approximation of a linear inlet, even if the fracture is not vertical and does not intersect the well on a very long distance. We might estimate this distance from the following argument. The flow rate observed at Soultz is about Q = 20 L/s. Thus, using a velocity of about u = 3.6 m/s and a fracture aperture equal to 3.6 mm implies that the well crosses such fractures over a cumulated length of about (neglecting the well radius):

$$L=\frac{Q}{uA}=\frac{20\times {10}^{-3}}{3.6\times 3.6\times {10}^{-3}}\approx 1.5\text{m}$$ |

_{x}/R

_{//}constant (here equal to 20.5): for instance, the results apply for l

_{x}= 690 m, A = 10 mm and ΔP/l

_{x}= −1.7 × 10

^{−4}bar/m (using the same fluid parameters), which leads to a Poiseuille velocity of about 0.46 m/s.

As we see, $\overline{T}$ is very inhomogeneous and also exhibits channeling. The chosen inlet temperature is T_{0} = 60 °C, the rock temperature is T_{r} = 200 °C (Fig. 7) and the fluid diffusivity is χ = 0.17 mm^{2}/s (corresponding to water at T = 100 °C, from Taine et Petit (2003). Note that the rock temperature will evolve over time in contrast to the one of our approximations. Indeed, the rock thermal diffusivity is about 1 mm^{2}/s which is larger than the fluid diffusivity (χ = 0.17 mm^{2}/s) but not sufficiently to be fully neglected.

However, $\overline{T}$ is rather different from a parallel plate solution. Indeed, the solution is not invariant in $\stackrel{\u02c6}{z}$. Different temperature profiles function of x are shown in Fig. 8. Two end-member types of behavior are plotted: temperature profiles at z = 960 m (curve iv) and z = 700 m (curve v). The temperature difference can be larger than 100 °C in the inlet region. Even rather far from the inlet, for example at x = 200 m, the temperature difference can still be of the order of 17 °C (200.0 °C for z = 960 m, and 183.4 °C for z = 700 m). The temperature field T(x,y,z = 700 m) for this set of parameters is shown in Fig. 3, where we see how the temperature evolves along the x-axis, away from the mean plane. Temperature profiles can be compared to the one obtained for a parallel plate model where plates are separated by the aperture A (curve iii) which reads from Eqs. (13) and (17): ${\overline{\overline{T}}}_{//}={\overline{T}}_{//}$.

Following this model (curve iii), the fluid should be at 199.7 ̊C at x = 200 m. If we compare ${\overline{\overline{T}}}_{//}$ to the averaged observed temperature $\overline{\overline{T}}$ (defined in Eq. (13)) (Fig. 8, curve i), we see that ${\overline{\overline{T}}}_{//}$ is not representative of the end-member types of behavior. We model $\overline{\overline{T}}$ by using an adapted parallel law ${\overline{\overline{T}}}_{\mathrm{mod}}$ (curve ii), which is an exponential law with a suitable thermal length R:

$${\overline{\overline{T}}}_{\mathrm{mod}}=({T}_{0}-{T}_{r})\mathrm{exp}\left(-\frac{x-{x}_{0}}{R}\right)+{T}_{r}$$ | (19) |

_{//}) and x

_{0}= −10 m. As we do a least square fit on a semi-log space to obtain the parameters R and x

_{0}, the beginning of the fit in the linear-linear space is not accurate. We see that the distance between wells (600 m in our case study) is about six times larger than R. However, owing to channeling the fluid temperature will not necessarily be in full equilibrium with the rock temperature at the out well. The thermal aperture is equal to Γ = 4.7 mm, which is rather different from the geometrical aperture A = 3.6 mm. A larger thermal aperture (compared to the geometrical one) means an inhibited thermalization on average.

### 3.2 Temperature estimation with few parameters

The knowledge of the spatial correlations rather than all the details of the geometrical aperture seems to be a key parameter to evaluate the hydraulic flow and the temperature of the fluid in a rough fracture. Indeed, the macroscopic geometrical aperture A brings too little information to characterize the heat exchange at the fracture scale. By contrast, it is impossible in particular for field measurements to know in detail the spatial variability of the local geometrical aperture A. Therefore, we propose to characterize the macroscopic geometrical properties with more than a single value, using several parameters describing the largest spatial variations. Numerically, it is possible to obtain them by filtering the fracture aperture field in the Fourier domain.

Fig. 9 shows the aperture field displayed in Fig. 4 once it has been filtered with the following criterion: only the Fourier coefficients fulfilling

$${\left(\frac{{k}_{x}}{2\pi}{l}_{x}\right)}^{2}+{\left(\frac{{k}_{z}}{2\pi}{l}_{z}\right)}^{2}<{2}^{2}$$ | (20) |

_{x}and k

_{z}are the wave vector coordinates along respectively the x and z-axes. Since the Fourier transform is discrete, it means that we only keep the Fourier components corresponding to the wave number $({n}_{x},{n}_{z})=(2\pi /({k}_{x}.{l}_{x}),2\pi /({k}_{z}.{l}_{z}))$ in {(0,0);(0,1);(1,0);(1,1)} (i.e. the average A and the first Fourier modes along x and z are left).

Let us assume that we only have these data available to evaluate the hydraulic flow and heat exchange. Using the same method and the same parameters as previously, we compute the pressure field corresponding to the filtered aperture field. In Fig. 10, we show the hydraulic aperture field we obtain. As we see, the high hydraulic aperture channel in the middle of the figure remains, while high frequency variations are removed. These large scale fluctuations of the hydraulic flow, computed from the knowledge of a very limited set of Fourier modes of the geometrical aperture, might be obtained from field measurements. Then, the corresponding temperature field shown in Fig. 11 is computed. The main features of the thermal field (Fig. 7) are still visible: the main channel is at the same position and the values are of the same order of magnitude. Despite small local differences, this substitution model gives a relevant description of what thermally happens.

## 4 Conclusions

We propose a numerical model to estimate the impact of the fracture roughness on the heat exchange at the fracture scale between a cold fluid and the hot surrounding rock. We assume the flow regime to be permanent and laminar. The numerical model is based on a lubrication approximation for the fluid flow (Reynolds equation). We also introduce a “thermal lubrication” approximation, which leads to a quartic profile of the temperature across the aperture. It is obtained by assuming that the in-plane convection is dominant with respect to the in-plane conduction (i.e. high in-plane Péclet number). The lubrication approximation implies also that the off-plane convection is neglected; subsequently, the heat conduction initiated by the temperature difference between the rock and the fluid is supposed to be the major off-plane phenomenon.

Our model shows that the roughness of the fracture can be responsible for fluid channeling inside the fracture. In this zone of high convection, the heat exchange is inhibited, i.e., the fluid needs a longer transport distance to reach the rock temperature. Spatial variability of the temperature is characterized on average by a thermal length and a thermal aperture.

In this article, we illustrate our modeling by a case study at the geothermal reservoir of Soultz-sous-Forêts, France, with a rough aperture which leads to inhibited thermal exchanges owing to a strong channeling effect, in the sense that the characteristic thermal length in this stationary situation is longer in rough fractures than in flat ones with the same transmissivity. Qualitatively, this can be attributed to the localization of the flow inside a rough fracture: most of the fluid flows through large aperture zones at faster velocity than the average, which leads to longer thermalization distances.

We performed simulations for about 1000 other aperture fields (not illustrated here) compatible with the macroscopic observations, of the same type as shown here. A general property holds for all these aperture fields: the thermal exchanges are always inhibited in rough fractures, compared to a fracture modeled by parallel plates with the same macroscopic hydraulic transmissivity (Neuville et al., submitted) (same H).

From the numerical solutions, we see that the mean geometrical aperture provides too little information to characterize the variability of the fluid flow and fluid temperature. In contrast, the knowledge of the dominant spatial variation of the geometrical aperture field (here obtained by keeping only the largest scale fluctuations using low pass Fourier filtering) provides interesting information about the spatial pattern of the hydraulic and thermal fields. The macroscopic spatial correlation of the aperture is shown to be an important parameter ruling the hydrothermal behavior. Note that we considered a self-affine model for the aperture roughness, but other types of geometrical descriptions of this roughness (given either by constraints from field measurements, or other kind of geometrical models), could be also considered using the type of simulations described here.

## Acknowledgements

We thank Albert Genter, Reinhard Jung Marion, Patrick Nami, Marion Schindler, Eirik G. Flekkøy, Stéphane Roux, Jose S. Andrade Jr. and Yves Méheust for fruitful discussions. This work has been supported by the EHDRA project, the REALISE program and the French Norwegian PICS project.

^{1} We compared our method to another algorithm based on a Lattice Boltzmann (LB) method, which does not reduce Navier-Stokes to a Stokes equation and does not hypothesize any lubrication approximation, in order to solve the velocity and temperature fields. The finite diffusivity of the rock is also taken into account. From those results, it appears that the analytical parabolic and quartic approximations (with the proper coefficients) of the respective fields are indeed consistent within a 5% error bar with the LB results.