# Crustal Fault Zones (CFZ) as Geothermal Power Systems: A Preliminary 3D THM Model Constrained by a Multidisciplinary Approach

Hugo Duwiquet, Laurent Guillou-Frottier, Laurent Arbaret, Mathieu Bellanger, Théophile Guillon, Michael Heap

► To cite this version:

Hugo Duwiquet, Laurent Guillou-Frottier, Laurent Arbaret, Mathieu Bellanger, Théophile Guillon, et al. Crustal Fault Zones (CFZ) as Geothermal Power Systems: A Preliminary 3D THM Model Constrained by a Multidisciplinary Approach. Geofluids, 2021, 2021, pp.8855632. 10.1155/2021/8855632. hal-03139470

HAL Id: hal-03139470

https://hal.science/hal-03139470

Submitted on 12 Feb 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research Article

# Crustal Fault Zones (CFZ) as Geothermal Power Systems: A Preliminary 3D THM Model Constrained by a Multidisciplinary Approach

Hugo Duwiquet<sup>1,2,3</sup> Laurent Guillou-Frottier,<sup>1,2</sup> Laurent Arbaret,<sup>2</sup> Mathieu Bellanger,<sup>3</sup> Théophile Guillon,<sup>1</sup> and Michael J. Heap<sup>4</sup>

Correspondence should be addressed to Hugo Duwiquet; hugo.duwiquet@cnrs-orleans.fr

Received 18 September 2020; Revised 7 December 2020; Accepted 16 December 2020; Published 4 February 2021

Academic Editor: Hang Lin

Copyright © 2021 Hugo Duwiquet et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

The Pontgibaud crustal fault zone (CFZ) in the French Massif Central provides an opportunity to evaluate the high-temperature geothermal potential of these naturally permeable zones. Previous 2D modeling of heat and mass transfer in a fault zone highlighted that a subvertical CFZ concentrates the highest temperature anomalies at shallow depths. By comparing the results of these large-scale 2D numerical models with field data, the depth of the 150°C isotherm was estimated to be at a depth of 2.5 km. However, these results did not consider 3D effects and interactions between fluids, deformation, and temperature. Here, field measurements are used to control the 3D geometry of the geological structures. New 2D (thin-section) and 3D (X-ray microtomography) observations point to a well-defined spatial propagation of fractures and voids, exhibiting the same fracture architecture at different scales (2.5 µm to 2 mm). Moreover, new measurements on porosity and permeability confirm that the highly fractured and altered samples are characterized by large permeability values, one of them reaching $10^{-12} m^2$. Based on a thermoporoelastic hypothesis, a preliminary 3D THM numerical model is presented. A first parametric study highlights the role of permeability, stress direction, and intensity on fluid flow. In particular, three different convective patterns have been identified (finger-like, blob-like, and double-like convective patterns). The results suggest that vertical deformation zones oriented at 30 and 70° with respect to the maximum horizontal stress direction would correspond to the potential target for high-temperature anomalies. Finally, a large-scale 3D numerical model of the Pontgibaud CFZ, based on THM coupling and the comparison with field data (temperature, heat flux, and electrical resistivity), allows us to explore the spatial geometry of the 150°C isotherm. Although simplified hypotheses have been used, 3D field data have been reproduced.

## 1. Introduction

As potential geothermal reservoirs, crustal fault zones (CFZ) are still largely unexplored and unexploited (e.g., South Erzgebirge deep fault zone, Germany [1]; Gravberg, Sweden [2, 3]; Badenweiler-Lenzkirch Sutur, Germany [4]; Dewey's Fault Zone, Canada [5]; Longmen Shan Fault Zone, China [6, 7]; Pasmajärvi Fault Zone, Finland [8]; and Luchlompal fault,

Russia [9]). The cost-effectiveness of a reservoir depends on permeability, heat, and fluid flow. These conditions are often related to a specific geodynamic context and the physical conditions of heat transport.

A CFZ is generally not a simple plane but rather has a complex geometry defined by a succession of fault and fracture intersections reflecting intense deformation [10–12]. CFZ are manifestations of strain localization and modify the

<footnote n="1"><sup>1</sup>BRGM, Georesources Division, Av. C. Guillemin, BP 36009, F-45060 Orléans Cedex 2, France</footnote>

<footnote n="2"><sup>2</sup>ISTO, UMR 7327, University of Orléans, CNRS, BRGM, 1A Rue de la Férollerie, 45071, Orléans, France</footnote>

<footnote n="3"><sup>3</sup>TLS-Geothermics, 91 Chemin de Gabardie, 31200 Toulouse, France</footnote>

<footnote n="4"><sup>4</sup>University of Strasbourg, CNRS, ENGEES, Institut Terre et Environnement de Strasbourg (ITES), UMR 7063, 5 Rue Descartes, Strasbourg F-67084, France</footnote>
mechanical properties of the crust above the Brittle-Ductile Transition (BDT). These heterogeneities, characterized by the presence of a succession of undeformed protolith, damage zones, and fault cores, will cause a significant variation in the permeability of the crust [13]. Permeability variation in the crustal domain has been widely discussed since the 1990s [14, 15]. Overall, the permeability of a CFZ is controlled by (i) 3D distribution of structural features such as brittle planes or lithological interfaces, (ii) the physical properties (e.g., Young's modulus and porosity) of lithologies along and adjacent to the fault zone, (iii) the state of regional stress, and (iv) more generally by the coupling between thermal (T), hydraulic (H), mechanical (M), and chemical (C) processes [10, 16–24]. At depths below the BDT, the deformation of granite leads to reductions in matrix porosity and, presumably, matrix permeability (e.g., [25]); however, these authors also note that microcracks can remain open in granite at high temperature (900°C) and high pressure (100 MPa), resulting in matrix permeabilities higher than that of the intact material.

Convective heat transfer is often invoked to explain temperature anomalies within faults [26–33]. Many parameters seem to condition the convective regime: the permeability ratio between the fault and its host, the CFZ permeability distribution, the fluid properties (viscosity and density variations), the CFZ thickness, and the depth and the dip of the CFZ [32–36]. The stress field is well known to influence CFZ permeability [20, 37–39] especially by influencing the void properties (fractures, pores, channels, and cracks). In fact, these voids under stress are expected to evolve over time and, according to several authors, favour shapes that will facilitate the circulation of fluids and stress propagation at different scales [40–45].

In order to estimate the geothermal potential of CFZ, it is necessary to consider a multidisciplinary approach [13, 15, 46–50]. Using such an approach (structural field measurements, laboratory porosity and permeability analysis, and numerical modeling) Duwiquet et al. [51] studied the geothermal potential of the Pontgibaud CFZ in the French Massif Central (FMC) and highlighted (i) the heterogeneous 2D distribution of permeability within the fault zone, (ii) CFZ with a near vertical dip as potential targets for high-temperature geothermal energy, and (iii) the position of the 150°C isotherm at a depth of 2.5 km. Although all these results were in agreement with previously acquired field data (geology, geophysics, and temperature measurements) [52, 53], the understanding of the 2D distribution of permeability has been made on a limited number of samples, and the modeling estimates have been made in 2D and based on TH coupling only. In addition, the structural framework poorly considered neighbouring faults to the Pontgibaud fault system. Using more field observations and more laboratory data, we aim to better characterize the geothermal system potential of the Pontgibaud CFZ in order to establish prospecting guides for geothermal exploration of CFZ worldwide.

Finally, to account for the mechanical aspect, we have performed a preliminary THM model where hypotheses are simple in order to focus on one single but important issue: could mechanical stresses affect the hydrothermal convective regimes? The objective of this simplified 3D THM model is

twofold: first, a simple synthetic model is dedicated to the demonstration of the effect of stress intensity and direction on fluid flow and thermal regime within the fault zone. Second, we implement mechanical effects within a 3D TH model of the Pontgibaud CFZ, with the unique objective of questioning consistency of the obtained results with field data. To our knowledge, THM models are generally applied to purely pressure-driven flow, where fluid density is constant and thus with no buoyancy-driven components (see, however, [54, 55]). However, Vallier et al. [24] used a similar approach as ours but with another objective being the study of mechanical effects on gravity anomalies, so that no information on the coupling between fluid flow and mechanical effects can be retrieved.

Thus, in this paper, we (i) provide additional qualitative and quantitative data necessary for the understanding of a hydrothermal system such as the Pontgibaud CFZ; (ii) insert these new data in a 3D geological model; (iii) explore by numerical modeling the effect of stress intensity, stress direction on fluid flow, and temperature anomalies; and (iv) evaluate in 3D the geothermal potential of the Pontgibaud CFZ.

## 1.1. Geological Setting

### 1.1.1. Geological Formation.
Bounded to the east and west by the Aigueperse Saint Sauves fault and the Sillon Houiller fault, the "La Sioule" licence (Figure 1) is largely affected by late Variscan tectonics. Metamorphic formations and the Gelles granite, which are intrusive, characterize the field. Magmatic and hydrothermal veins are widely represented.

Covering a large part of the mapped region, metamorphic formations (green in Figure 1) can be separated into different lithologies. Biotite and sillimanite paragneiss extends southwards near the Miouze river and northeastwards to the Villelongue region. Orthogneiss associated with biotite and sillimanite paragneiss are present sporadically. Cordierite gneisses are more widely represented, but it is difficult to assign a cartographic boundary. Rare amphibolitic levels are interspersed in the cordierite gneisses. As an intrusive in metamorphic series, the porphyroid granite of Gelles (red in Figure 1) extends southwest of Pontgibaud, between the Sioule Valley to the East and the village of Tortebesse to the West. To the North, the last granitic points appear immediately to the west of La Goutelle.

The presence of thermomineral springs (such as Pranal and Ceyssat) is an indicator of present-day fluid circulation (Figure 1). Temperature measurements have been carried out in the Pontgibaud borehole between a depth of 55 and 120 m (International Heat Flow Commission database, http://ihfc-iugg.org/; borehole number FR198). From these measurements, a geothermal gradient of $41°C/km^{-1}$ has been inferred, and a heat flow of $110 mW·m^{-2}$ was estimated. Magnetotelluric campaigns have been conducted in the Pontgibaud fault zone area by TLS-Geothermics and IMAGIR between 2015 and 2017 [50, 51]. Two anomalies located in the brittle crust at depths of 3 and 8 km could be linked to the presence of magma, clays (smectites),
<figure><img src="image_2.png" /></figure>

FIGURE 1: Geological and geothermal features in the “La Sioule” licence area (black contour). FMC: French Massif Central. Heat flow and geothermal gradient (red crosses [56, 57], International Heat Flow Commission). Heat production [58] and resistivity profile (black dashed line [51]). P: Pontgibaud; LG: La Goutelle; V: Villelongue; T: Tortebesse. Coordinate system: WGS84 Pseudo-Mercator EPSG:3857 (after the 1/50,000th geological map published by BRGM).

and/or hydrothermal fluids. These two anomalies appear to
be related to the Pontgibaud fault zone.

1.1.2. *Structural Setting.* In order to understand the spatial organization of potential drains of some geological structure, it is necessary to first consider the tectonic evolution of the study area (Figure 2). We will then investigate the effects of the direction of the principal stress on the fluids flow in preexisting fractures.

(1) Famennian to late Namurian (365 Ma-315 Ma)

From Famennian to the early Visean, the maximum main horizontal stress axis is rotated 40° clockwise from NNW-SSE to NNW-SSW. During the Visean, the maximum main stress axis spatially varied from the vertical (extension, detachment of Argentat) to the horizontal (strike-slip, sinistral kinematic of the “Sillon Houiller”). Then, during the late Namurian, the maximum main horizontal stress shortening axis
<figure><img src="image_3.png" /><figcaption>FIGURE 2: Paleostress orientation in the “La Sioule” licence since the Tournaisian (360 Ma). Blue arrows indicate the extensional direction. Red arrows indicate the compressional direction (modified after [59, 65, 66]).</figcaption></figure>

continued to rotate clockwise until it reached a NE-SW direction (Figure 2).

### (2) Early to late Stephanian (305-295 Ma)

During the lower and middle Stephanian, N-S compression created a system of conjugated N140°-170°E sinistral and N10°-50°E dextral with N80°-120°E reverse faults. This deformation episode was responsible for the NE-SW sinistral fault of the “Sillon Houiller.” During the transition from the middle to the upper Stephanian, the axis of rotation of the compression changed from N-S to NW-SE (Figure 2). This rotation activated or reactivated dextral conjugated faults N 120°E and sinistral N160-180°E. During the late Stephanian, the compression became E-W and resulted in a strong folding of coal deposits along the Sillon Houiller fault [59].

### (3) Permian to Jurassic (295-205 Ma)

Following Stephanian compression, the N-S/NE-SW Permo-Triassic extension is represented in the basement by normal faults, tension cracks, and veins [59–62]. On the Cévennes border, this N-S extension followed NW-SE directions inherited from the Stephanian, forming NE faults such as the Cévennes fault [59, 63]. At the Trias-Lias border, a major metallogenic event occurred in the Massif Central, with the formation of a large number of veins and

stratiform deposits, whose age was determined at about 194 Ma [64].

### (4) Jurassic to Cretaceous (205-96 Ma)

During the transition from the Triassic to the Jurassic, a NE-SW extension occurred in the Rodez and in the western Cévennes, creating normal synsedimentary faults [59]. The beginning of the NW-SE extension was accompanied by significant metallogenic and thermal phenomena. Studies show that the appearance of N-S- to E-W-oriented barite and fluorite veins in Morvan, Quercy, and Lozère can be attributed to the Early Jurassic [61, 66, 67]. The Jurassic-Cretaceous period corresponds mainly to the development of the Liguro-Provençal Ocean and then the Atlantic Ocean. The beginning of the Jurassic is characterized by a radical change in the directions of extension, changing from a period marked by a strong late extension to an extensive period where the stresses are much lower, causing small deformations. This N-S oriented extension episode was responsible for the appearance of normal N80°-130°-oriented faults. The rapid rotation of the direction of the minimum $σ_3$ and intermediate $σ_2$ stresses around the maximum subvertical stress $σ_1$ led to a new N-S-oriented paleogeography that persisted in the early and mid-Jurassic (Figure 2). At the beginning of the Late Jurassic, the Dogger
structure was destroyed by another rotation of the minimum $\sigma_3$ and intermediate $\sigma_2$ stresses around the maximum vertical stress $\sigma_1$. The previous E-W extension disrupted the organized N-S paleogeography elevations of the Middle Jurassic.

### (5) Cretaceous to late Pliocene (96-2.5 Ma)

The establishment of a rift characterizes the Cenozoic in Europe: the European Cenozoic Rift System (ECRIS). This large rift system (about 1000 km) is associated with a series of transect faults throughout the Variscan massif. This produced a radical change in the stress regime throughout the Massif Central (Figure 2). After a long period of Permian-Mesozoic extension, a N-S compressive tectonic regime was established, the effects of which were visible from England to Bohemia [68–71].

In the Massif Central and its borders, the Late Eocene-Oligocene extension generated basins, grabens, and half-grabens bounded by normal faults, in which several hundred to thousand metres of continental sediments were deposited. The faults limiting these basins are generally N-S oriented but some NE-SW and NW-SE trends also exist. While the extension affected the entire Massif Central and its border, the deformation was concentrated near the regional fractures inherited from the Variscan orogeny.

During the Oligocene, the Western European plate was affected by a system of N-S-trending grabens forming the Western European Rift. This continental rift can be divided into several sections linked by transforming regional faults with a NE direction [71, 72]. Among the N-S grabens, the Limagne half-graben, 15 km east of the PFZ, hosts a temperature anomaly that reaches 160°C at a depth of 3 km [73].

### (6) Present-day stress regime

Several seismological stress-state evaluations were performed around the studied area. Mazabraud et al. [74] have determined five focal mechanisms from regional records below the Pulvérières and St. Gervais granites, located 20 km north of the PFZ. They exhibit reverse, normal, and strike-slip mechanisms. The maximum horizontal stress axis is oriented close to NW-SE. More recently, local stress states were evaluated during the exploration of the geothermal licence based on the seismological records acquired between 2015 and 2016 by TLS-Geothermics using 11 triaxial broadband seismometers CMG40 (see [53] for more details about this geophysical acquisition). Data were processed in 2018 by MAGNITUDE (a Baker Hughes company) to (i) identify microearthquake events, (ii) perform a focal mechanism analysis, and (iii) evaluate the stress state. 25 microearthquake events were identified and localized (mainly to the north of the studied area), and 14 focal mechanisms were determined (between -2 and -12 km asl). Confidence is very good for six

Focal Mechanism Solutions (FMS) (good amplitudes fit at each station and focal spheres are well constrained), is good for three FMS, and is more uncertain for five FMS. They exhibit reverse (three FMS), normal (two FMS), and strike-slip regimes (nine FMS). The strike-slip regime occurred for all FMS with a high degree of confidence, for 2/3 of FMS with a good degree of confidence, and for 1/5 of FMS with a low degree of confidence. The nature of the stress tensor for each FMS quality set (see [74]) was then estimated by assuming that the following hypotheses are true: (i) the stress field is uniform within the studied area and (ii) the tangential traction on the fault plane is parallel to the slip direction. In each of the three cases, they obtain a strike-slip regime (i.e., with a subvertical intermediary stress axis and subhorizontal maximum and minimum stress axes). Directions of the maximum stress axis are bracketed between N135° and N150° for the high quality set and between N120° and N135° for the other two cases (Figure 2).

## 2. Materials and Methods

The multidisciplinary approach used here corresponds to a method widely used in oil and uranium exploration, CO₂ capture and storage studies, and reservoir engineering. It combines a field approach (structural measurements and sampling of the Peyrouses 1 borehole), a laboratory approach (2D and 3D observations and measurements of permeability and connected porosity), and, finally, numerical modelling. The goal of our study is to integrate these field and laboratory measurements and observations into a 3D numerical model that is based on a Thermal, Hydraulic, and Mechanical (THM) coupling method using Comsol Multiphysics™.

### 2.1. *Structural Observations.*
The structural measurements, performed in the field, make it possible to better characterize the intensity of the deformation in and around the Pontgibaud fault zone and also to put the geometry and architecture flow pathways into a 3D perspective. Sampling of the Peyrouses 1 borehole located within the Pontgibaud fault system provides a better understanding of the spatial distribution of pores and permeable fractures in this hydrothermal system, through thin-section observations and X-ray microtomography.

### 2.2. *Laboratory Observations and Measurements.*
The samples used for laboratory measurements and observations were recovered from the Peyrouses 1 borehole (Figure 3), located within the PFZ (see Figure 1). These samples represent different facies: altered, fractured, altered and fractured, and intact facies (anatexite and microgranite). Thin-section observations on samples with different degrees of fracturing and alteration and at different scales help to characterize the petrology of the Pontgibaud hydrothermal system. 3D observations using a Phoenix nanotom microtomograph have been used at the Institut des Sciences de la Terre d'Orléans (ISTO). X-ray microtomography is a high-resolution, nondestructive 3D imaging technique. The samples were placed between a tungsten filament and a detector. During
<figure><img src="image_4.png" /><figcaption>FIGURE 3: Distribution of lithologies and location of the various samples taken from the Peyrouses 1 borehole. Mineralization consists of mispickel, pyrite, tourmaline, pyromorphite, and galena. For the missing samples (in black on the figure), the description of the drilling done by BRGM in 1981 was used.</figcaption></figure>

acquisition, the sample rotated 360° and the tungsten filament emitted X-rays onto the sample (using a voltage and current of 110 kV and 85 µA, respectively). The detector measured the attenuation of the emitted signal. A series of 2D images were produced, with pixel values dependent on the intensity of the emitted radiation and on the attenuation recorded by the detector. This attenuation is a function of the density and atomic number of the material through which the rays pass [76]. A 3D reconstruction of this object has been

done, and each voxel reflects the average attenuation of the traversed material. However, the density contrast between two mineral phases may bias the information. A Beam-Hardening Correction (BHC) limits the loss of this information. After the reconstruction, the 3D volume was processed with the Volume Graphics VG Studio Max visualization software. This imaging technique provides high-resolution images (2.5 to 10 µm) and allows us to understand how fluids can flow in these heterogeneous materials. X-ray microtomography is a common method used to characterize oil [77], uranium [78], and geothermal reservoirs [79].

In order to quantify the heterogeneity of the system, permeability measurements have been performed at the Institut de Physique du Globe de Strasbourg (IPGS) in Strasbourg (France). For this purpose, 11 cylindrical samples (with a diameter of 25 mm and length between 25 and 40 mm) were prepared from an intact sample, an altered sample, a fractured sample, and an altered and fractured sample. The crystal size of these samples was small when compared to their diameter, as is recommended for samples used for measurements of permeability [79]. These samples were first oven dried under a vacuum at 40°C. Their connected porosities were calculated using the connected volume measured by a helium pycnometer (Micromeritics AccuPyc II 1340) and the bulk sample volume measured using digital calipers. Gas permeability was measured using a benchtop nitrogen parameter (see the schematic diagram in [81, 82]). All measurements were performed under a confining pressure of 1 MPa (samples were first left for 1 h at this pressure for microstructural equilibration) using the steady-state method. Volumetric flow rates were measured (using a Bronkhorst gas flowmeter) for six pressure differentials (measured using a Keller pressure transducer) to calculate permeability using Darcy's law and to check whether the Forchheimer and/or Klinkenberg corrections are required (for more details on this protocol, see [83, 84]). The 3D numerical modeling discussed below will use these laboratory data.

## 2.3. *Numerical Approach*

Considering a very simple case, we will first try to understand the effects of stress direction and stress intensity and the effects of permeability on fluid flow, intensity, and depth of temperature anomalies. Then, we perform a simple preliminary study on the natural case of the Pontgibaud CFZ in order to get consistent results with field data. We first describe below the governing equations together with physical parameters. Then, the idealized simple model is described before the 3D geological model.

The poroelasticity interface of Comsol Multiphysics™ combines the regular fluid flow formulation with linear elastic solid mechanics of the porous medium matrix. Poroelasticity coupling means that the deformation of the porous medium affects fluid flow, while changes in volumetric stress will in turn affect the amount of motion and heat transport. This interface includes an expression of the stress tensor as a function of the strain tensor and the Biot-Willis coefficient [85].

Mixing mass conservation of the fluid and momentum balance yields that the poroelastic model is governed by the following relation:
$$ \rho_f S \frac{\partial p_f}{\partial t} + \nabla \cdot (\rho_f \vec{v}) = Q_m - \rho_f \alpha_B \frac{\partial}{\partial t} \varepsilon_{vol}, \quad (1) $$

where $\rho_f$ is the fluid density (kg·m⁻³), S is the poroelastic storage coefficients (Pa⁻¹, see Equation (7)), $p_f$ is the fluid pressure (Pa), t is time, $Q_m$ is the fluid source term (kg·m⁻³·s⁻¹), $\alpha_B$ is the Biot-Willis coefficient (-), and $\varepsilon_{vol}$ is the volumetric strain of the porous matrix (-); Darcy's law sets that $\vec{v}$ is equal to

$$ \vec{v} = -\frac{k}{\mu_f} (\nabla p_f - \rho_f \vec{g}), \quad (2) $$

where k is the permeability of the medium (m²), $\mu_f$ is the fluid dynamic viscosity (Pa·s), and g is acceleration due to gravity (m·s⁻²).

The governing equation for the momentum balance model is

$$ -\nabla \cdot \sigma = \rho g. \quad (3) $$

Here, $\sigma$ (Pa) is the stress tensor and $\rho$ (kg·m⁻³) is the total density. Notice that negative compression convention is considered.

The components of the strain tensor depend on the displacement vector $\vec{u}$ (m). Under the small strain assumption, we have

$$ \varepsilon = \frac{\nabla \vec{u} + T \nabla \vec{u}}{2}. \quad (4) $$

The stress-strain relationship for linear materials relates the stress tensor and the strain tensor $\varepsilon$ through Hooke's law. Accounting for the hydromechanical coupling, we have

$$ \sigma = 2\mu\varepsilon + 3\lambda\varepsilon_V\delta - \alpha_B p_f \delta, \quad (5) $$

for isotropic materials. Here, $\delta$ is the second order identity tensor (-), $\varepsilon_V = Tr(\varepsilon)/3$ is the volumetric strain (-), and $\mu$ and $\lambda$ are Lamé's parameters (Pa) and are a function of Young's modulus E (Pa) and Poisson's ratio $\nu$ (-) by the following relationship:

$$ \lambda = \frac{E\nu}{(1 + \nu)(1 - 2\nu)}, \quad (6) \\ \mu = \frac{E}{2(1 + \nu)}. $$

The storage coefficient S (in Equation (1)) depends on the porosity $\varepsilon_p$ (-), Biot coefficient $\alpha_B$ (-), fluid compressibility $\chi_f$ fixed at $4.8 \times 10^{-10}$ Pa⁻¹ [86], and matrix compressibility $\chi_m$ (Pa⁻¹):

$$ S = \varepsilon_p \chi_f + \chi_m (\alpha_B - \varepsilon_p) (1 - \alpha_B). \quad (7) $$

The matrix compressibility is fixed with Young's modulus and Poisson's coefficient by the following relationship:

$$ \chi_m = \frac{E/3}{1 - (2 \times \nu)}. \quad (8) $$

Consequences of Equations (7) and (8) on numerical values of S are the following: the coefficient is time-dependent through porosity variations and the obtained values vary between $10^{-11}$ and $10^{-10}$ Pa⁻¹. These values are consistent with the compressibility values provided by Freeze and Cherry [87].

The poroelastic material model obeys Equation (3) which describes changes in the stress tensor $\sigma$ and porous matrix deformation due to boundary conditions and changes in pore pressure. Equation (3) defines a state of static equilibrium because the changes in the solid equilibrate quickly, implying that there are no time-dependent terms. Still the time rate of change in strain $(\partial/\partial t)\varepsilon_{vol}$ appears as a coupling term in Equation (1) because the solid equation becomes transient when solved simultaneously with a time-dependent flow model.

Details of the coupling equations between Darcy's law and Fourier's law can be found in Guillou-Frottier et al. [88]. For our 3D experiments, realistic fluid and rock properties have been used. The fluid density law corresponds to fitting experimental data for pure water and for temperature and pressure ranges of 0-1000°C and 0-500 MPa, respectively (NIST Chemistry WebBook [89, 90]):

$$ \begin{aligned} \rho_L(T, P) = & 1006 + (7.424 \times 10^{-7} \times P) + (-0.3922 \times T) \\ & + (-4.441 \times 10^{-15} \times P^2) + (4.547 \times 10^{-9} \times P \times T) \\ & + (-0.003774 \times T^2) + (1.451 \times 10^{-23} \times P^3) \\ & + (-1.793 \times 10^{-17} \times P^2 \times T) + (7.485 \times 10^{-12} \times P \times T^2) \\ & + (2.955 \times 10^{-6} \times T^3) + (-1.463 \times 10^{-32} \times P^4) \\ & + (1.361 \times 10^{-26} \times P^3 \times T) + (4.018 \times 10^{-21} \times P^2 \times T^2) \\ & + 7.372 \times 10^{-15} \times P \times T^3) + (5.698 \times 10^{-11} \times T^4), \end{aligned} \quad (9) $$

where T is the temperature (°C) and P is the pressure (Pa). The dynamic viscosity of the fluid is a function of temperature T, such that [91, 92]

$$ \mu_f(T) = 2.414 \times 10^{-5} \cdot \exp \left( \frac{570}{T + 133} \right). \quad (10) $$

Some rock and fluid properties are set as constant for the parametric study and in the large-scale numerical modeling of the Pontgibaud system. Their values are summarized in Tables 1, 2, and 3 (see Supplementary Material).

Permeability is a critical parameter in the functioning of geothermal systems and will be our main adjustment variable. We will distinguish between the permeability of the basement and the permeability of the fault zone, whose variations will be differentiated. In the basement, and up to the BDT, the permeability will be subject to the lithostatic pressure and will therefore be considered as decreasing with depth. It will be characterized by the following equation [93]:
<figure><img src="image_5.png" /><figcaption>FIGURE 4: Geometry, meshing, and boundary conditions of the parametric study. Numerical set-up is detailed in the text.</figcaption></figure>

$$K_b(z) = K_{b_0} \exp \left( \frac{z - 800}{\delta} \right). \qquad (11)$$

Here, $K_b(z)$ (m²) is the depth-dependent permeability of the basement and $K_{b_0}$ (m²) is the maximum permeability (present at the surface) taken here at $10^{-16}$ m². $z$ (m) is the depth below sea level, and $\delta$ (m) defines the decrease in permeability with depth, taken here at 2500 m, following the models of [51, 88, 93].

When considering the fault zone permeability $K_f$, we have to consider the depth variation (as Equation (9)) but we must also consider the lateral variations observed and measured in the field and in the laboratory. This function of permeability is similar to that already used in [51], where the permeability varies in 2D along the $z$ and $x$ axes with a decrease in permeability with depth. In order to consider the permeability in 3D, we add a “$y$” coordinate to this function and we then obtain

$$K_f(x, y, z) = K_{f_0} \times \left[ \exp \frac{z - 800}{\delta} \right] \times \left[ 101 + 100 \times \sin \left( 2 \times \pi \times \frac{(x + y + z)}{\lambda} \right) \right], \qquad (12)$$

where $K_f(x, y, z)$ is the permeability of the fault dependent on the three space directions and $K_{f_0} \times 201$ will be the maximum permeability at the surface. The choice of 100 and 101 means that the permeability will vary over two orders of magnitude. The permeability here will vary according to a sinusoid applied to the three directions $x, y, z$. The term $\lambda$ corresponds to the wavelength of the sinusoid. In order to approximate the observed and measured permeability variation (see below), we take a small value of $\lambda$.

The numerical simulations have been performed with Comsol Multiphysics™ in two different stages, illustrated in Figures 4 and 5, respectively. The first one is aimed at studying the effects of mechanics on fluid dynamics, through simplified and idealized hypotheses, while in the second one, previous geological observations and measurements are incorporated in a large-scale 3D THM model.

The first stage consists of investigating the effects of the fault direction, of the fault permeability, and of the stress value on the fluid flow, on the intensity of the temperature anomaly and on its depth. To do this, a series of numerical experiments have been performed with three variables: directions of the vertical deformation zone that vary from 0° to 90° (with a 10° step), stress values that vary from 0 to 100 MPa (with a 10 MPa step), and permeability values that vary from $1 \times 10^{-16}$ to $4 \times 10^{-14}$ m² (with a $10^{-16}$ m² step) (Figure 4). The mesh consists of 175,789 tetrahedra, with a maximum size of 200 m for low permeability values and a maximum size of 20 m for high permeability values. The stress is applied on the red side (see Figure 4); the other three lateral borders have roller conditions (no perpendicular displacement). At the base of the model, an embedding condition (no displacement) and a heat flow of 80 mW·m⁻² are applied. The vertical sides are thermally and hydraulically insulated, and at the surface, a temperature and pressure of 10°C and $10^5$ Pa are imposed, respectively. A positive pressure evolution with depth is imposed as an initial condition.

The initial thermal conditions of these models are defined from the stationary state obtained for a pure conductive regime. At $t=0$, values of permeability are assigned to the different compartments and the fluid flow is initiated.

The second stage concerns the numerical modeling of the Pontgibaud hydrothermal system in 3D (Figure 5). Here, the 3D geometry is constrained by the field observations and measurements. Indeed, based on the tectonic history (see Figure 2) and the measurements acquired in the field, different fault families have been individualized. The complete mesh consists
<figure><img src="image_6.png" /><figcaption>FIGURE 5: Geometry, meshing, and boundary conditions for numerical modelling of the Pontgibaud crustal fault zone. Numerical set-up is detailed in the text. Stress profile of the Chassoles borehole (45.38°N, 3.09°E) (from [95]) integrated in Comsol Multiphysics™.</figcaption></figure>

of 244,475 tetrahedra, with a maximum size of 200 m and a maximum size of 10 m in most permeable zones. The mesh is refined for the highest permeable zones as well as for fault and lithology intersections. The result of this numerical modelling will be compared with field data, such as heat flux and geothermal gradient, as well as resistivity data from [53]. Considering the thermal relaxation time of magmatic chambers [94], calculations were performed up to a time of $t_0 + 15$ kyrs.

One of the critical points in setting up regionally scaled mechanical models are the boundary conditions [96, 97] that result from the far-field stresses. The overall stress orientation is often hinted from the integration of several measurements in a given region (such as well tests and GPS measurement). In our case, the present-day stress orientation in the surroundings of the Pontgibaud CFZ is approximately N145° (see section “structural setting,” above). As far as
magnitudes are concerned, few or no indications usually
exist, and our case is no exception. In most cases, numerical
models include the boundary magnitudes as unknowns
which are estimated through sensitivity analyses [98] or
inversion studies [96]. In our case, we propose to take
advantage of *in situ* measurements from the Chassoles bore-
hole located at 60 km south-east of Pontgibaud ([95]
Figure 4). This well was drilled for the “Géologie Profonde
de la France” research program. The stress field has been
determined by hydraulic testing. The well encountered only
orthogneiss and intersected several subvertical fault zones at
depths of 180 and 630 m. A total of 22 tests were run, nine
between 340 and 470 m and 13 between 734 and 853 m.
Many tests yielded fractures. At Chassoles, down to 7 km,
the maximum stress is vertical but nearly equal to $\sigma_H$. This
stress profile has been implemented as boundary conditions
for the 3D model of the Pontgibaud CFZ. The heterogene-
ities in rock properties and structural features affect the
stress distributions [99], and the extrapolation of measure-
ments to other locations and/or depths appears to be a
strong assumption. Still, regarding the high number of
unknowns and parameters in our model, we decided to keep
these observations as inputs in our first computations for
the sake of simplicity.

**3. Results**

3.1. *Structural Characterization.* The north-south structures (Figure 6, A, C, D, and H) represent fairly clear topographic boundaries and underline the axis of Pb-Zn mineralized veins. The dips of these structures vary between 60 and 90°. The alignment of this axis is often interrupted by transverse structures and masked by the abundance of lavas in the eastern part. Although masked by the alluvium of the Sioule Valley, these north-south structures (Figure 6, C and D) underline a homogeneous alignment over 40 km.

E-W structures also mark the topography (Figure 6, E and K). With a dip between 60 and 90°, field observations and the geological map show an occurrence of these structures every 10 km in the southern part of the FMC and near the Sillon Houiller.

The NW-SE (Figure 6, G, I, and J) structures are quite frequent and often overlap the E-W faults without inducing significant offset.

The NE-SW (Figure 6, B, D, and F) structures can be seen
widely on the geological map, from the west of the Sillon
Houiller to Limagne. They are particularly well expressed
on the Gelles granite. Mapping observations show that these
faults are not too affected by N-S faults. The dip of these fault
planes are mostly vertical.

To sum up, the Pontgibaud CFZ is a large deformation
zone (>3 km) extending over 40 km, mainly characterized
by north-south structures, and crossed by E-W, NE-SW,
NW-SE, and some other NNE-SSW structures. These field
observations have been used to determine the 3D geometry
of our Pontgibaud large-scale numerical model. Before that,
we will visualize and quantify the permeability distribution
of the PFZ.

3.2. *Thin-Section Observations (2D).* The following section summarizes the main lithologies present as well as the void space arrangement, for different scales of 2D observations. In these images, the voids are represented in white and the fractures in red. In order to make sure that the voids present are not caused by thin-section manufacture, we have highlighted the voids where secondary mineral phases have precipitated.

The facies of the samples are essentially fractured and
altered (Figures 7(b), 7(c), 7(d), and 7(f)). Quartz, feldspars,
biotites, and muscovites are mainly found. There are also
secondary mineral phases (such as tourmaline, pyrite, arse-
nopyrite, scorodite, and sericite) and accessory minerals
(such as apatite and zircon). The least altered and unfrac-
tured sample is shown in Figure 7(a) (LPA view). Quartz
are most often xenomorphic and can form clusters of small
crystals or develop over a wide range with an automorphic
tendency, giving the rock a porphyritic texture. Biotite is also
observed and generally automorphic and shows an onset of
chloritization. We also observe the polysynthetic macle of
plagioclase in Figure 7(a) (LPA view).

Within these altered and fractured facies, voids are pres-
ent within cores of quartz mineralization (Figure 7(b)) but
also at the walls of these mineralizations (Figures 7(b) and
7(c)). While some fractures have voids through which fluids
can flow (Figure 6, D and E), other fractures are completely
filled with clays, secondary mineral phases, and/or oxides.
In addition, fractures can also be cemented by cataclasites
(Figure 7(f)). These cataclasites consist of angular quartz
clasts, plagioclase, sericite clays, and some calcites. Due to
their similar petrological properties, most of these angular
clasts should originate from the basement.

Most voids and secondary mineral phases appear to be
associated with the development of fractures and cracks.
Their interconnection could serve as potential drains for fluid
circulation. The architecture of this interconnection at differ-
ent scales could have a key role on the ability of the system to
allow fluids to circulate.

Figures 8(a)–8(d) show images of the same sample but at
progressively smaller scales (from 2 mm to 50 µm). These
zoomed images are focused at the core of a fracture network
in order to understand the spatial distribution of voids in the
central part of the fracture network. Figures 8(e)–8(h) again
represent the same sample but at progressively smaller scales
(from 500 µm to 50 µm). These zoomed images are focused at
the extremities of a fracture network in order to see how frac-
tures and pores spread to the extremities of these networks.

At a scale of 2 mm, Figure 8(a) shows macroscopic voids.
They are distributed within the quartz vein and also at the
walls of this mineralization. The secondary minerals and
other smaller voids appear to form a network of fractures
(red line) that extend away from the quartz vein, forming a
dendritic network. At a scale of 500 µm (Figure 8(b)), other
voids can be seen, still marked by secondary mineralization,
which also follows the dendritic fracture network. By increas-
ing the zoom again, we observe voids (Figure 8(c)) that we
could not distinguish in Figure 8(b). These voids are still
present within this dendritic network and are still marked
by secondary mineralization. By increasing the zoom again,
<figure><img src="image_7.png" /><figcaption>FIGURE 6: Simplified structural map, showing the main fractures present in the “La Sioule” licence and the regional stress orientation (Schmidt’s projection on the lower hemisphere. Data from TLS-Geothermics [95]). Coordinate system: WGS84 Pseudo-Mercator EPSG:3857 (after the 1/50,000th geological map published by BRGM).</figcaption></figure>

up to 50 µm, we still notice voids that we could not distinguish at the previous scale. These voids also follow a dendritic pattern. We find an interconnection of fractures at different scales, which form a dendritic fracture network. The precipitation of secondary mineral phases inside this dendritic fracture network shows that fluids have circulated.

Figure 8(e) shows another sample with a fracture network organization identical to Figure 8(a). Fractures marked by secondary mineralization, or by the presence of voids, appear to propagate spatially. By zooming in on Figure 8(e), we find voids that we could not distinguish previously (Figure 8(f)). These voids seem to follow the fractures and are also marked

by the presence of late mineralization. Zooming in again on Figure 8(h), we now see a distribution of fractures that was not visible in Figure 8(g) and with an architecture not yet observed. This propagation of fractures across the box is in the form of rectangular drainage patterns (Figure 8(h)).

These flow patterns are widely present in natural drainage systems. These observations give us an overview of the spatial representation of the fractures and voids. We will now see if these structures are connected in 3D.

3.3. *X-Ray Microtomography Observations (3D)*. The visualization of fractures and 3D voids in a sample taken from the
<figure><img src="image_8.png" /></figure>

FIGURE 7: Facies, textures, and composition details of the Pontgibaud hydrothermal system. (a) LPA thin-section view of least altered and unfractured samples. Others (b–e) show altered and fractured facies which are much more present in the samples. (c) is a zoom of a part of (b) and shows voids present at smaller scales. (f) Crosspolar thin-section view of cataclasite. Abbreviations from Whitney and Evans, [100]: Pl: plagioclase; Bt: biotite; Qz: quartz; Ser: sericite. Other abbreviations: Cat: cataclasite; Ox: oxyde. The thin section come from samples retrieved from the Peyrouses 1 borehole, located within the Pontgibaud CFZ (see location in Figure 1).

Pontgibaud CFZ (see location in Figure 1) are shown in Figures 9(a)–9(d). The sample is characterized by a fractured part as well as a matrix part. Samples were cored at the edge between the fractures and the matrix (Figure 9(b)) as well as in the matrix (Figures 9(c) and 9(d)).

The fractures are shown in red (Figure 9(b)), and the voids are shown in yellow (Figures 9(b)–9(d)). The fractures extend over the entire length of the sample and are heterogeneously distributed (Figures 9(b)–9(d)) in space. With a resolution of 10 µm, this 3D illustration of a complex architecture of the fracture network shows an interconnection of the voids present in these fractures (Figure 9(b)) that could facilitate the circulation of fluids. The 3D visualization of the matrix (Figure 9(c)) highlights the voids with a resolution of 10 µm. At this resolution, one observes globular clusters that appear to be interconnected. With a finer resolution of 2.5 µm, we can see more pores present in the matrix (Figure 9(d)). This

resolution allows us to confirm the presence of voids that form a planar structure.

Thin-section (2D) and X-ray microtomography (3D) observations allow us to observe the architecture of the fracture networks and the distribution of pores at different scales both within the fracture networks and in the matrix. The porous flow structures of the Pontgibaud hydrothermal system have multiple scales and a nonuniform distribution of the space available for fluid circulation. The anisotropic nature of the distribution of permeability in the different dimensions of space has thus been quantified by measurements of permeabilities and connected porosities.

### 3.4. Porosity and Permeability Measurements.
The data show that there is a general trend of increasing permeability as a function of increasing porosity (Figure 10). Connected porosity values range from 6 to 22% (Figure 10). In most
<figure><img src="image_9.png" /><figcaption>FIGURE 8: Fracture design for two samples ((a–d) and (e–h), respectively) and for smaller observation scales. Fractures underlined by red lines are marked by secondary mineral phases (mainly tourmaline and oxyde) and voids.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 9: 3D reconstruction of the matrix and fractures of a sample from the Peyrouses 1 borehole. Fractures are shown in red and voids in yellow. The resolution of these images is 10, 10, and 2.5 µm ((b-d), respectively).</figcaption></figure>

cases, the altered samples (with and without fractures) have the highest connected porosity values (around 20%), and the intact samples have the lowest porosities (between 6 and 8%). Permeability values range from $2 \times 10^{-18}$ to $1 \times 10^{-12}$ m$^2$. Altered samples with fractures have the highest permeability values (between $10^{-13}$ and $10^{-12}$ m$^2$). Samples that are either altered or contain a fracture are characterized by intermediate permeabilities, between $10^{-16}$ and $10^{-14}$ m$^2$. The intact samples have the lowest permeability anisotropy, which is most pronounced in the samples that contain fractures (with and without alteration) (Figure 10). For example, the unaltered sample containing a fracture (Figure 10) is a factor of 20 more permeable in the X direction (blue dot in Figure 10) than in the Y direction (green dot in Figure 10) (despite their similar porosities). The permeability of the altered sample containing a fracture (Figure 10) is three orders of magnitude more permeable in the Y direction (green cross in Figure 10) than in the X direction (blue cross in Figure 10), despite their similar porosities.

### 3.5. Numerical Approach.
More than 2500 numerical experiments were run to understand the effects of fault direction, stress intensity, and permeability on the intensity and depth of the positive temperature anomaly. Figure 11 synthesizes

these results. This figure presents the value of the positive temperature anomalies, as well as their depths, as a function of the stress intensity (horizontal axis), the angle with respect to the stress direction (vertical axis), and the permeability ratio between the fault and the host rocks defined as

$$R = \frac{K_{f \max}}{K_{b_0}}. \qquad (13)$$

This permeability ratio compares the maximum value of the permeability of the fault ($K_{f \max}$) with that of the basement ($K_{b_0}$) (here fixed at $10^{-16}$ m$^2$). Other ($K_{b_0}$) values were tested, but higher values do not correspond to the field data, and lower values of permeability do not change the results.

### 3.5.1. Numerical Parametric Study.
By considering all the directions of the deformation zone and all the stress values, we can see that the largest temperature anomalies and the smallest depths will be concentrated for values of R [50-200] and R [200-400] (Figure 11(a)). By varying the direction of the deformation zone with respect to the stress, for permeability ratios R [50-400], the higher temperature anomalies and the smallest depths will be concentrated in directions
<figure><img src="image_11.png" /></figure>

FIGURE 10: Permeability as a function of connected porosity for altered, fractured, altered and fractured, and unaltered and unfractured samples from the Peyrouses 1 borehole (see location in Figure 1). \*Data from Duwiquet et al. [51]. The X, Y, and Z directions (respectively, in blue, green, and red) correspond to the sampling direction for each core for the same facies. Measurement uncertainty corresponds to the size of the symbol.

of 30 and 70°. This tendency is, however, less visible for R [1-50]. By varying the stress, we notice that for R [1-50], the intensity of the stress does not significantly influence the value of the positive temperature anomaly as well as its depth. However, trends seem to emerge for stresses of 30 and 90 MPa. These stress values seem to concentrate the highest temperature anomalies for the directions of 30 and 70° relative to the stress as well as for R [200-400]. For values of R [50-200], we do not observe significant effects of stress on the depth and intensity of the temperature anomaly.

Finally, for R [50-400], the deformation zone directions of 30 and 70° relative to stress will concentrate the highest temperature anomalies at the shallowest depths. For these directions, stresses of 30-90 MPa will have the highest temperature anomalies.

Figures 11(b) and 11(c) show the fluid flow for R [50-200]. For Figure 11(b), the direction of the deformation zone is oriented at 30° to the stress, which is 30 MPa. The fluid flow (marked by the black arrows) is defined by an upward movement in the center and downward movements at both ends. The bulging of the 70°C isotherm (shown in red in Figure 11(b)) is the consequence of this fluid circulation, which follows a “blob-like” pattern. For Figure 11(c), the fault direction is 70° with respect to the stress, which is 50 MPa. Fluid flow is also defined by an upward movement in the center and downward movements at both ends. However, compared to Figure 11(b), the downward movements are more concentrated around the upward movement. The morphology of the 70°C isotherm (shown in red in Figure 11(c)) is the consequence of

this fluid circulation which follows a “finger-like” pattern. As can be seen, there seems to be different types of convection patterns for different stress and fault direction values. Note that the irregular morphology of the red surfaces may be due either to the 3D heterogeneity of the permeability or to the effects of stress direction and/or intensity. To solve for this important issue, convective patterns have been systematically observed for all experiments.

Observations of the convective figures from the numerical simulations make it possible to distinguish three zones as a function of the stress intensity, fault direction, and the value of the temperature anomaly (Figure 12). The brown zone and the green zone correspond to a ratio R [1-200]. The red zone corresponds to a ratio R [200-400]. For this permeability ratio, the temperature anomalies will be between 88 and 150°C. For the lowest permeability ratio (R [1-200]), temperature anomalies are between 7 and 88°C. For these values, two convective regimes are observable. Firstly, for stress values below 40 MPa and for a fault direction below 50°, the blob-like pattern described in Figure 11 is found. This convective regime occurs along two dimensions of the space (along the *z* and *x* axes). For stresses greater than 40 MPa and for a fault direction greater than 50°, the finger-like pattern described in Figure 11 is preferred. This convective regime is expressed along two dimensions of the space (along the *z* and *x* axes). For the highest permeability ratios (R [200-400]), we find a continuous upward and downward motion in all three dimensions of space. Circulation of warmest fluids is along the *z* and *y* axes, while the
<figure><img src="image_12.png" /><figcaption>Example of convective regimes</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c)</figcaption></figure>

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

FIGURE 11: (a) Positive temperature anomalies and their depth for different permeability ratios (black, green, and red) as a function of fault direction (vertical axis) and stress intensity (horizontal axis). The underlined values show the trends of temperature anomalies and their depths for fault direction of 30 and 70° and for different stress values. (b, c) show the fluid flow for the same R ratio (green), but for different fault direction and stress. Red surface illustrates the morphology of the 70°C isotherm. (b) shows a "blob-like" pattern. (c) shows a "finger-like" pattern.

circulation of the coldest fluid is along the *z* and *x* axes. This circulation of fluids occurs according to a "double-like" pattern.

By considering a simple geometry and idealized hypotheses, it can be seen that the stress direction, the stress intensity, and the permeability ratio between the fault and the basement are parameters that have an influence on fluid flow and consequently on the intensity and depth of the temperature anomalies. Considering now a more complex geometry, constrained by field data, we can try to see if similar trends as

those found previously (finger-like, blob-like, and double-like convective patterns) can be retrieved in a 3D THM numerical model.

### 3.5.2. Large-Scale Numerical Modeling Validation.
The results of the numerical modeling show that for a maximum permeability imposed on the Pontgibaud deformation zone of $1 \times 10^{-14} \text{ m}^2$, the geothermal gradient would be 41.5°C/km (Figure 13(a)) and the heat flux would be $110 \text{ mW} \cdot \text{m}^{-2}$ (Figure 13(b)), at the location (Figure 1) where the measured
<figure><img src="image_15.png" /><figcaption>FIGURE 12: Diagram of convective regimes. The distorted planes represent the isotherms. The arrows represent the circulation of fluids. In the brown area are the blob-like convection patterns. In the green zone are the finger-like convection patterns. In the red zone are the double-convective patterns.</figcaption></figure>

and estimated field data are 41°C/km and 110 mW·m⁻². The resistivity profiles show a low resistivity anomaly at a depth of 1 km (Figure 13(b)) and a larger low resistivity anomaly at a depth of 3-4 km. If the origin of this resistivity may be discussed, the comparison between the resistivity profiles and the temperature anomalies shows that the temperature anomaly within the Pontgibaud deformation zone is between 2 and 3 km deep. In the middle part of the Pontgibaud deformation zone and at the intersections of the different fault families (N-S, NE-SW, NW-SE, and E-W), the 150°C isotherm (in red on Figure 13(a)) would occur at a depth of 2.3 km.

Permeability variation (Figure 13(c)) imposed on large-scale numerical simulations corresponds to observations made in the field, in 2D thin-section images, and in 3D X-ray microtomographic images.

## 4. Discussion

4.1. *Origin and Effect of Heterogeneities.* Thin-section (2D) and X-ray microtomography (3D) observations have revealed spatially propagating fracture networks in which fluids can flow through pores at different scales, from 2.5 µm to 2 mm (Figures 7, 8, and 9). This fracture propagation has been

studied on mineralized systems (St. Salvy, France, and Jebel Aouam, Morocco) and in which tectonic comminution processes have been identified [101–103]. Microscopic scale pores have also been observed in the walls of some mineralization. The most frequent process in the brittle part of the crust is hydrofracturing. This process depends on the temporal variation of the fluid pressure [104–106], causing a “break up” which may be at the origin of some voids. This process has been studied on sites with F (El Hammam-Morocco), Ba (Dreislar, Germany), and Ag and Au (Creede, Colorado, USA) mineralization [102]. Fluid circulation can also cause the dissolution of mineral phases in the basement. Dissolution may be accompanied by a loss of volume and thus cause the formation of pores within the matrix [90]. While the consequences of these processes can be observed in 2D and 3D, other processes can cause the opening or closing of the voids (such as volume reduction, volume expansion, and pore collapse), thus influencing the capacity of the system to allow fluids to circulate (Figure 14).

Measurements of permeabilities and connected porosities from samples from Pontgibaud show ranges of values between 2 × 10⁻¹⁸ and 1 × 10⁻¹² m² and 6 and 22%, respectively (Figure 10). These values are within known ranges of
<figure><img src="image_16.png" /><figcaption>FIGURE 13: Pontgibaud large-scale numerical simulation results (a) for a maximum permeability imposed on the Pontgibaud CFZ of $1 \times 10^{-14} \text{ m}^2$ (c) are compared with the resistivity profiles, the geothermal gradient, and the heat flux (b). Fluid velocities vary from $10^{-8}$ to $10^{-10} \text{ m} \cdot \text{s}^{-1}$. This permeability variation is the result of Equation (10).</figcaption></figure>

permeability for rocks in a basement context [1]. Our data also highlight the heterogeneities of the permeability of a hydrothermal system, which could be caused by the processes of tectonic comminution, hydrofracturing, and mineral phase dissolution, processes that can form structures and pores at different scales. The permeability of a rock mass will be very much controlled by the fracture density and, importantly for hydrothermal systems, the proportion of open fractures [82, 84, 107, 108]. Our permeability data also show that alteration can increase permeability (Figure 10). Previous studies on volcanic hydrothermal systems have also shown that hydrothermal alteration can greatly influence matrix permeability [109–113]. For these natural hydrothermal systems, permeability should not be approximated to a single value but should rather be represented as pipes or channels embedded in the neighbouring system [114].

According to Bejan and Lorente [115], architectures at different scales and arranged in a hierarchical pattern would promote fluid flow. These natural drainage systems (Figures 8 and 9) may have lower flow resistance than a similarly porous volume. The constructal theory [40] allows us to describe the natural tendency of flow systems to generate multiscale shapes and structures favouring fluid flow. It has also been shown that the transition between heat transfer at the nanoscale and more

conventional scales is governed by the smallest dimensions and by their spatial heterogeneities [116]. The origin of these heterogeneities could be the stress (with tectonic comminution) and the water (with fluid pressure). Finally, with knowledge of current parameters, the geometry of the Pontgibaud CFZ might correspond to an optimal one, which facilitates efficiency of convective flow.

## 4.2. Effect of Stress Direction and Stress Intensity on the Fluid Flow.
Multiple parameters can influence fluid flow, and the 3D behaviour of convective cells can be complex [117]. López and Smith [31] showed that temperature anomalies can be directly influenced by the permeability ratio between the fault and the host. 3D modeling based on TH coupling has shown the influence of this permeability ratio on the morphology of the temperature anomalies [33]. Our parametric study (Figures 11 and 12) made it possible to approach three factors influencing fluid flow: stress intensity, stress direction, and the permeability ratio between the fault and the basement. For high permeability ratios between $R$ [200-400], the poroelastic effects due to stress state variation have no significant effect on fluid flow pattern. However, for lower permeability ratios $R$ [1-200], changes in the fluid flow pattern can be observed depending on the stress direction and stress intensity.
<figure><img src="image_17.png" /><figcaption>FIGURE 14: (a) Three-dimensional conceptual model of the Pontgibaud CFZ as a geothermal reservoir potential. Details of the voids present at different scales and in which fluids can flow. (b) The processes that may be at the origin of these voids (tectonic comminution, fluid-assisted brecciation, and basement alteration) are illustrated. (c) Comparison of the large-scale modeling results with field data shows that fluids could flow through a blob-like convective pattern and locate the 150°C isotherm at 2.3 km depth.</figcaption></figure>

### 4.3. Positive Temperature Anomaly on the Pontgibaud CFZ.
As shown in Figure 13, the results of Pontgibaud large-scale numerical modeling compared with heat flux, geothermal gradient, and resistivity data show that the 150°C isotherm could be located at a depth of 2.3 km within the Pontgibaud CFZ. This upwelling of the 150°C isotherm is the result of fluid circulation through a blob-like convective pattern. This will produce a temperature anomaly within the Pontgibaud hydrothermal system, but lateral thermal diffusion is observed within the fault zone itself, thus defining a larger geothermal reservoir.

Indeed, considering a simplified 3D geometry, this lateral temperature diffusion has been highlighted in recent numerical models of fluid flow around faulted zones [33].

The Pontgibaud CFZ is located at 45° to the main regional stress (Figure 14). In this direction, the convective regime diagram (Figure 12) predicts blob-like convective pattern fluid flows for R ratios [1-200], as well as for stress intensities below 40 MPa. In view of these parameters, the result of the convective regime diagram seems to be consistent with the validation of the large-scale numerical model of the Pontgibaud CFZ in which blob-like convective pattern fluid flows would be present.

If the consistency between general concepts from the very simplified numerical modeling is effectively validated by natural system results, then it would be possible to use the convective pattern diagram in order to characterize convection regimes of a natural system according to a tectonic regime.

Figures 13(a) and 14(c) show a positive temperature anomaly that extends widely on both sides of the Pontgibaud CFZ. The origin of this anomaly would be a blob-like convective pattern. Both blob-like and finger-like convective patterns are found for low permeability ratios (see Figure 12). Finger-like

convective patterns are characterized by isotherms covering a smaller space than blob-like convective patterns. Consequently, the volume of the temperature anomaly will be less in a finger-like convective pattern than in a blob-like one. As we have seen, the parameters influencing the formation of these both convective pattern are the intensity and the direction of the stress. Thus, in a geothermal exploration context, in geological settings where the permeability is quite low, it would be advantageous to target areas where the intensity of the stress is quite low and where the direction of the stress is less than 50° with respect to the deformation zone. In this case, there will be more widespread temperature anomalies.

For high permeability zones, the intensity and direction of the stress are no longer parameters influencing the convective regime under thermoporoelastic conditions. The double-convective pattern gives higher temperature anomalies than the two patterns previously described (ΔT > 87°C). Although the temperature anomaly is more intense, the deformation of the isotherms is more restricted than for blob-like convective patterns. Thus, in areas of high permeability, the temperature anomaly will be the most intense, but the volume covered will be less than for blob-like convective patterns.

### 4.4. Limitation and Perspectives.
The theory of thermoporoelasticity considers the two effects of temperature change and fluid pressure on the elastic behaviour of rocks. This theory could be further treated in mechanical terms, in which slip and dilation tendency could be calculated [118, 119]. In fact, some fault segments subjected to critical stresses would act as conduits for the circulation of fluids [120]. As a consequence, some preferential flow paths can appear under the action of mechanics.

The synthetic model illustrates the impact of mechanics with one boundary stress only. The effects are expected to
be greatly emphasized under anisotropic boundary conditions because of the shearing conditions they introduce.

Permeability is widely recognized as a central parameter in the understanding of a geothermal system. It is also a factor that evolves in space and time. Stresses, fluid overpressure, thermal cracking, mineral dissolution, and precipitation are just a handful of examples that can affect the permeability of a geothermal system [121–125]. All of these phenomena help to explain how a system that is initially impermeable evolves into a geothermal reservoir potential over time. Another method to improve such 3D models would be to be able to consider permeability as a dynamic parameter that varies in response to different geological processes, including varying tectonic regime as illustrated by the Pontgibaud CFZ ().

The parametric studies were realized by considering the faults as vertical, faults which from the 2D parametric study on TH coupling showed that the highest temperature anomalies concentrate at the lowest depth. However, considering a THM coupling, it would be interesting to be able to consider the dip of the structures. It would then be possible to highlight stress directions, permeability ratios, stress intensity, and dips as favourable targets for high temperature geothermal energy. In general, geothermal exploration must consider the intensity of the stress, the direction of the stress, and the variability of permeability.

Finally, all these results are only applicable for large deformation zones in which temperature anomalies are induced by buoyancy-driven flow and not in zones where there is significant topography. Moreover, we only consider one fluid phase (liquid water). If we were to consider two phases, liquid and gas, the poroelastic effects on the compressibility of the gas could be an important process [126]. If the compressibility of the fluid can also influence the elastoplastic behaviour [127], the thermal dependence of compressibility might change the thermal distribution of the system [128].

## 5. Conclusions

Based on the results of this multidisciplinary study, we speculate that crustal fault zones represent potential reservoirs for high-temperature geothermal energy. Field measurements collected at the Pontgibaud crustal fault zone highlight a large N-S deformation zone, with a with a dip of 60° to the east and intersected by different families of NE-SW, NW-SE, and E-W faults. Laboratory observations show drainage structures with mainly dendritic geometry. The origin of this geometry here could be the fluid pressure and/or tectonic commutation. Such flow structures are widely present in different natural systems and have been suggested by the structural theory as optimal geometries for fluid flow. All these observations and measurements have been inserted in 3D numerical models with THM coupling. The first simplified numerical study allowed us to conclude that

(i) for the simplified used boundary conditions, vertical deformation zones, having an angle of 30° and 70° with respect to the boundary stress, will concentrate the highest temperature anomalies, at the lowest depths

(ii) different convective patterns could exist: finger-like convective patterns, blob-like convective patterns, and double-like convective patterns. These three convective regimes depend on parameters such as the intensity and direction of the stress, the permeability ratio between the basement, and the deformation zone

In a second step, large-scale numerical simulations of the Pontgibaud crustal fault zone were performed. Recalling that no calibration was conducted, these computations were ran essentially as a demonstration of their ability to study natural phenomena while incorporating field observations. Interestingly, they already show consistent results with field data. These results suggest that the 150°C isotherm is at a depth of 2.3 km and could therefore represent an interesting target for high-temperature geothermal exploration. This tendency will be more thoroughly studied in future works by using more constrained models (with, e.g., sensitivity analyses). Also, these results need to be confirmed by drilling by the operator. If the drilling confirms these results, then it could be used on other study zones such as Badenweiler-Lenzkirch Sutur (Germany), Longmen Shan Fault Zone (China), and Pasmajärvi Fault Zone (Finland).

## Data Availability

Sample position, structural measurements, permeability and connected porosity measurements are included in the paper (Figures 3, 6, and 10). The physical properties of fluids and rocks are detailed in Tables 1, 2, and 3 (see Supplementary Material).

## Conflicts of Interest

The authors declare that there are no conflicts of interest for the publication of this article.

## Acknowledgments

We sincerely thank the association "La Route des Mines Dômes-Combrailles" as well as the "Musée de la Mine" for allowing the sampling of the Peyrouses 1 borehole. We sincerely thank Laurent Bailly for his help on the thin-section observations and for his constructive remarks. We warmly thank Prof. Fabien Magri for these valuable suggestions and remarks. We are grateful to Séverine Caritg, Eric Fournier, Anne-Sophie Serrand, Antoine Armandine Les Landes, Simon Lopez, Arnold Blaisonneau, Bastien Hermant, and Jean-Michel Ars, for their fruitful discussions and relevant advice. The financial support was provided by TLS-Geothermics, ISTO, and BRGM and by the "GERESFAULT" project (ANR-19-CE05-0043).

## Supplementary Materials

The supplementary materials contain three tables. They group together the physical and thermal properties of rocks and fluids inserted in the numerical modeling. (Supplementary Materials)
## References

[1] P. Achtziger-Zupančič, S. Loew, and A. Hiller, "Factors controlling the permeability distribution in fault vein zones surrounding granitic intrusions (Ore Mountains/Germany)," *Journal of Geophysical Research: Solid Earth*, vol. 122, no. 3, pp. 1876-1899, 2017.

[2] C. Juhlin and H. Sandstedt, *Storage of Nuclear Waste in Very Deep Boreholes: Feasibility Study and Assessment of Economic Potential. Pt. 1 and 2*, Swedish Nuclear Fuel and Waste Management Co., 1989.

[3] B. Lund and M. D. Zoback, "Orientation and magnitude of in situ stress to 6.5 km depth in the Baltic Shield," *International Journal of Rock Mechanics and Mining Sciences*, vol. 36, no. 2, pp. 169-190, 1999.

[4] O. Brockamp, A. Schlegel, and K. Wemmer, "Complex hydrothermal alteration and illite K-Ar ages in Upper Visean molasse sediments and magmatic rocks of the Variscan Badenweiler-Lenzkirch suture zone, Black Forest, Germany," *International Journal of Earth Sciences*, vol. 104, no. 3, pp. 683-702, 2015.

[5] C. Bieber, D. Chorley, W. Zawadzki, and J. Reinson, "Hydrogeologic data collection and development of conceptual models to predict mine inflow quantity and quality at Diavik diamond mine, NWT," in *Proceedings of the Sea to Sky Geotechnique 2006 Conference*, Vancouver, BC, Canada, 2006.

[6] T. X. Yang, Q. Duan, J. Chen, and M. J. Dekkers, "Rock magnetic expression of fluid infiltration in the Yingxiu-Beichuan fault (Longmen Shan thrust belt, China)," *Geochemistry, Geophysics, Geosystems*, vol. 17, no. 3, pp. 1065-1085, 2016.

[7] Q. Duan, X. Yang, and J. Chen, "Hydraulic properties of a low permeable rupture zone on the Yingxiu-Beichuan fault activated during the Wenchuan earthquake, China: implications for fluid conduction, fault sealing, and dynamic weakening mechanisms," *Tectonophysics*, vol. 721, pp. 123-142, 2017.

[8] A. E. Ojala, J. Mattila, M. Markovaara-Koivisto, T. Ruskeenemi, J. P. Palmu, and R. Sutinen, "Distribution and morphology of landslides in northern Finland: an analysis of postglacial seismic activity," *Geomorphology*, vol. 326, pp. 190-201, 2019.

[9] Y. A. Kozlovsky, *The Superdeep Well of the Kola Peninsula. Exploration of the Deep Continental Crust*, Springer-Verlag, New York, 1987.

[10] F. M. Chester and J. M. Logan, "Implications for mechanical properties of brittle faults from observations of the Punchbowl fault zone, California," *Pure and Applied Geophysics*, vol. 124, no. 1-2, pp. 79-106, 1986.

[11] S. E. Schulz and J. P. Evans, "Spatial variability in microscopic deformation and composition of the Punchbowl fault, southern California: implications for mechanisms, fluid-rock interaction, and fault morphology," *Tectonophysics*, vol. 295, no. 1-2, pp. 223-244, 1998.

[12] C. G. Sammis, A. J. Rosakis, and H. Bhat, *Effects of off-fault damage on earthquake rupture propagation: experimental studies*, Mechanics, Structure and Evolution of Fault Zones, Birkhäuser Basel, 2009.

[13] D. R. Faulkner, C. A. L. Jackson, R. J. Lunn et al., "A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones," *Journal of Structural Geology*, vol. 32, no. 11, pp. 1557-1575, 2010.

[14] S. J. Martel, "Formation of compound strike-slip fault zones, Mount Abbot quadrangle, California," *Journal of Structural Geology*, vol. 12, no. 7, pp. 869-882, 1990.

[15] V. F. Bense, T. Gleeson, S. E. Loveless, O. Bour, and J. Scibek, "Fault zone hydrogeology," *Earth-Science Reviews*, vol. 127, pp. 171-192, 2013.

[16] R. H. Sibson, "Brecciation processes in fault zones: inferences from earthquake rupturing," *Pure and Applied Geophysics*, vol. 124, no. 1-2, pp. 159-175, 1986.

[17] W. T. Parry, P. N. Wilson, and R. L. Bruhn, "Pore-fluid chemistry and chemical reactions on the Wasatch normal fault, Utah," *Geochimica et Cosmochimica Acta*, vol. 52, no. 8, pp. 2053-2063, 1988.

[18] J. D. Byerlee, "Model for episodic flow of high-pressure water in fault zones before earthquakes," *Geology*, vol. 21, no. 4, pp. 303-306, 1993.

[19] R. L. Bruhn, W. T. Parry, W. A. Yonkee, and T. Thompson, "Fracturing and hydrothermal alteration in normal fault zones," *Pure and Applied Geophysics*, vol. 142, no. 3-4, pp. 609-644, 1994.

[20] C. A. Barton, M. D. Zoback, and D. Moos, "Fluid flow along potentially active faults in crystalline rock," *Geology*, vol. 23, no. 8, pp. 683-686, 1995.

[21] S. A. Miller, A. Amos, and D. L. Olgaard, "Earthquakes as a coupled shear stress-high pore pressure dynamical system," *Geophysical Research Letters*, vol. 23, no. 2, pp. 197-200, 1996.

[22] C. H. Scholz, *The Debate on the Strength of Crustal Fault Zones*, The Mechanisms of Earthquake Faulting, 2002.

[23] A. A. L. Landes, T. Guillon, M. Peter-Borie, A. Blaisonneau, X. Rachez, and S. Gentier, "Locating geothermal resources: insights from 3D stress and flow models at the upper Rhine Graben scale," *Geofluids*, vol. 2019, 24 pages, 2019.

[24] B. Vallier, V. Magnenet, J. Schmittbuhl, and C. Fond, "HM modeling of gravity anomalies related to deep hydrothermal circulation at Soultz-sous-Forêts (France)," *Geothermal Energy*, vol. 8, no. 1, pp. 1-21, 2020.

[25] M. Violay, M. J. Heap, M. Acosta, and C. Madonna, "Porosity evolution at the brittle-ductile transition in the continental crust: implications for deep hydro-geothermal circulation," *Scientific Reports*, vol. 7, no. 1, pp. 1-10, 2017.

[26] R. P. Lowell, "Circulation in fractures, hot springs, and convective heat transport on mid-ocean ridge crests," *Geophysical Journal International*, vol. 40, no. 3, pp. 351-365, 1975.

[27] F. Magri, S. Möller, and N. Inbar, "2D and 3D coexisting modes of thermal convection in fractured hydrothermal systems - Implications for transboundary flow in the Lower Yarmouk Gorge," *Marine and Petroleum Geology*, vol. 78, pp. 750-758, 2016.

[28] J. Faulds, M. Coolbaugh, V. Bouchot, I. Moeck, and K. Oguz, *Characterizing structural controls of geothermal reservoirs in the Great Basin, USA, and Western Turkey: developing successful exploration strategies in extended terranes*, 2010.

[29] H. D. Murphy, "Convective instabilities in vertical fractures and faults," *Journal of Geophysical Research: Solid Earth*, vol. 84, no. B11, pp. 6121-6130, 1979.

[30] B. Vallier, V. Magnenet, J. Schmittbuhl, and C. Fond, "THM modeling of hydrothermal circulation at Rittershoffen geothermal site, France," *Geothermal Energy*, vol. 6, no. 1, p. 22, 2018.

[31] D. L. López and L. Smith, "Fluid flow in fault zones: analysis of the interplay of convective circulation and topographically driven groundwater flow," *Water Resources Research*, vol. 31, no. 6, pp. 1489-1503, 1995.
[32] C. Forster and L. Smith, "The influence of groundwater flow on thermal regimes in mountainous terrain: a model study," *Journal of Geophysical Research: Solid Earth*, vol. 94, no. B7, pp. 9439-9451, 1989.

[33] L. Guillou-Frottier, H. Duwiquet, G. Launay, A. Taillefer, V. Roche, and G. Link, "On the morphology and amplitude of 2D and 3D thermal anomalies induced by buoyancy-driven flow within and around fault zones," *Solid Earth*, vol. 11, no. 4, pp. 1571-1595, 2020.

[34] V. I. Malkovsky and A. A. Pek, "Timescales for reaching steady-state fluid flow in systems perturbed by formation of highly permeable faults," *Geofluids*, vol. 4, no. 4, pp. 253-258, 2004.

[35] V. I. Malkovsky and F. Magri, "Thermal convection of temperature-dependent viscous fluids within three-dimensional faulted geothermal systems: estimation from linear and numerical analyses," *Water Resources Research*, vol. 52, no. 4, pp. 2855-2867, 2016.

[36] P. Achtziger-Zupančič, S. Loew, A. Hiller, and G. Mariethoz, "3D fluid flow in fault zones of crystalline basement rocks (Poehla-Tellerhaeuser ore field, Ore Mountains, Germany)," *Geofluids*, vol. 16, no. 4, pp. 688-710, 2016.

[37] R. H. Sibson, "Structural permeability of fluid-driven fault-fracture meshes," *Journal of Structural Geology*, vol. 18, no. 8, pp. 1031-1042, 1996.

[38] T. Ito and M. D. Zoback, "Fracture permeability and in situ stress to 7 km depth in the KTB scientific drillhole," *Geophysical Research Letters*, vol. 27, no. 7, pp. 1045-1048, 2000.

[39] J. V. Rowland and R. H. Sibson, "Structural controls on hydrothermal flow in a segmented rift system, Taupo Volcanic Zone, New Zealand," *Geofluids*, vol. 4, no. 4, pp. 259-283, 2004.

[40] A. Bejan and S. Lorente, "Constructal theory of generation of configuration in nature and engineering," *Journal of applied physics*, vol. 100, no. 4, p. 5, 2006.

[41] A. Bejan, *Shape and Structure, from Engineering to Nature*, Cambridge university press, 2000.

[42] A. Bejan, "Constructal-theory network of conducting paths for cooling a heat generating volume," *International Journal of Heat and Mass Transfer*, vol. 40, no. 4, pp. 799-816, 1997.

[43] A. H. Reis, *Constructal Theory: From Engineering to Physics, and How Flow Systems Develop Shape and Structure*, 2006.

[44] A. F. Miguel, "Natural flow systems: acquiring their construc-tal morphology," *International Journal of Design & Nature and Ecodynamics*, vol. 5, no. 3, pp. 230-241, 2010.

[45] A. Bejan and J. P. Zane, "Design in nature," *Mechanical Engineering*, vol. 134, no. 6, pp. 42-47, 2012.

[46] M. Oda, "An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses," *Water Resources Research*, vol. 22, no. 13, pp. 1845-1856, 1986.

[47] J. S. Caine, J. P. Evans, and C. B. Forster, "Fault zone architecture and permeability structure," *Geology*, vol. 24, no. 11, pp. 1025-1028, 1996.

[48] V. F. Bense and M. A. Person, "Faults as conduit-barrier systems to fluid flow in siliciclastic sedimentary aquifers," *Water Resources Research*, vol. 42, no. 5, 2006.

[49] L. Micarelli, A. Benedicto, and C. A. J. Wibberley, "Structural evolution and permeability of normal fault zones in highly porous carbonate rocks," *Journal of Structural Geology*, vol. 28, no. 7, pp. 1214-1227, 2006.

[50] J. Alcalde, I. Marzán, E. Saura et al., "3D geological characterization of the Hontomin CO₂ storage site, Spain: multidisciplinary approach from seismic, well-log and regional data," *Tectonophysics*, vol. 627, pp. 6-25, 2014.

[51] H. Duwiquet, L. Arbaret, L. Guillou-Frottier, M. J. Heap, and M. Bellanger, "On the geothermal potential of crustal fault zones: a case study from the Pontgibaud area (French Massif Central, France)," *Geothermal Energy*, vol. 7, no. 1, p. 33, 2019.

[52] M. Bellanger, "High temperature geothermal resources of crustal fault zones: a dedicated approach," in *79th EAGE conference and exhibition*, Paris, France, 2017.

[53] J. M. Ars, P. Tarits, S. Hautot, M. Bellanger, O. Coutant, and M. Maia, "Joint inversion of gravity and surface wave data constrained by magnetotelluric: application to deep geothermal exploration of crustal fault zone in felsic basement," *Geothermics*, vol. 80, pp. 56-68, 2019.

[54] Y. Wang, T. Li, Y. Chen, and G. Ma, "A three-dimensional thermo-hydro-mechanical coupled model for enhanced geothermal systems (EGS) embedded with discrete fracture networks," *Computer Methods in Applied Mechanics and Engineering*, vol. 356, pp. 465-489, 2019.

[55] Y. Shi, X. Song, J. Li, G. Wang, R. Zheng, and F. YuLong, "Numerical investigation on heat extraction performance of a multilateral-well enhanced geothermal system with a discrete fracture network," *Fuel*, vol. 244, pp. 207-226, 2019.

[56] F. Lucazeau and G. Vasseur, "Heat flow density data from France and surrounding margins," *Tectonophysics*, vol. 164, no. 2-4, pp. 251-258, 1989.

[57] G. Vasseur, G. Robert, B. Feuga, and G. Bienfait, "Groundwater flow and heat flow in an area of mineral springs," *Geothermics*, vol. 20, no. 3, pp. 99-117, 1991.

[58] F. Lucazeau, *Flux de chaleur, production de chaleur et évolution géodynamique récente du Massif Central Français*, Thèse de doctorat, 1981.

[59] J. L. Blès, D. Bonyjoly, C. Castain, and Y. Gros, "Successive post-Variscan stress fields in the French Massif Central and its borders (Western European plate): comparison with geodynamic data," *Tectonophysics*, vol. 169, no. 1-3, pp. 79-111, 1989.

[60] G. Santouil, *Tectonique et microtectonique de la distension permienne et de l'évolution post-triasique dans les bassins de Lodève, Saint-Afrique et Rodez (France Sud-Ouest)*, Thèse Doct. 3e cycle. U.S.T.L, Montpellier, 1980.

[61] D. Bonijoly, *Étude structurale et minéralisations d'une plate-forme carbonatée, le Quercy*, Université d'Orleans, 1981.

[62] J. L. Blès and J. L. Bles, *Apport de l'analyse structurale à la connaissance des gîtes filoniens*, 1982.

[63] D. Bonijoly and H. Germain, "Histoire tectonique post-hercynienne du bassin d'Alès (Gard): chronologie des déformations et contrôle structural des minéralisations barytiques," *Editions du Bureau de recherches géologiques et minières, Service géologique national*, vol. 76, 1984.

[64] J. C. Baubron, *Nouvelles datations K/Ar sur des filons à quartz et fluorine du Massif Central français*, 1980.

[65] M. Bellanger, B. Hermant, S. Galibert, and J. L. Auxiètre, "Fault-controlled hydrothermal system associated with major Crustal Fault Zone: future drilling target to assess the deep geothermal potential-The Sioule license project," in *Proceedings European Geothermal Congress Den Haag*, The Netherlands, 2019.
[66] M. Jébrak, *Le filon des Farges et les minéralisations à barytine, fluorine et galène du district d'Ussel dans leur cadre géologique (Massif Central Francais)*, Thèse Doct. 3e cycle, University Orléans, 1978.

[67] C. O. Valette, *Karsts et filons & fluorine dans le faisceau synclinal du Morvan: le gisement d'Argentolle (Saône-et-Loire, France)*, Thèse Doct. 3e cycle. University Orléans-Dot. B. R. G. M., 1983.

[68] P. A. Ziegler, "Late Cretaceous and Cenozoic intra-plate compressional deformations in the Alpine foreland-a geodynamic model," *Tectonophysics*, vol. 137, no. 1-4, pp. 389-420, 1987.

[69] M. Mattauer and J. L. Mercier, *Microtectonique et grande tectonique*, vol. 10, pp. 141-161, 1980.

[70] J. Letouzey and P. Trémolières, *Paleo-stress fields around the Mediterranean since the Mesozoic derived from microtectonics: comparison with plate tectonic data*, vol. 115, pp. 261-272, 1980.

[71] F. Bergerat, *Déformations cassantes et champs de contraintes tertiaires dans la plate-forme européenne*, p. 315, 1985.

[72] J. P. Gélard, "A fracturation de la Bourgogne méridionale, essai d'interprétation et implications tectoniques régionales," *Revue de géographie physique et de géologie dynamique Paris*, vol. 20, no. 1, pp. 5-27, 1978.

[73] P. Calcagno, C. Baujard, L. Guillou-Frottier, A. Dagalier, and A. Genter, "Estimation of the deep geothermal potential within the Tertiary Limagne basin (French Massif Central): an integrated 3D geological and thermal approach," *Geothermics*, vol. 51, pp. 496-508, 2014.

[74] Y. Mazabraud, N. Béthoux, J. Guilbert, and O. Bellier, "Evidence for short-scale stress field variations within intraplate central-western France," *Geophysical Journal International*, vol. 160, no. 1, pp. 161-178, 2005.

[75] J. W. Gephart and D. W. Forsyth, "An improved method for determining the regional stress tensor using earthquake focal mechanism data: application to the San Fernando earthquake sequence," *Journal of Geophysical Research: Solid Earth*, vol. 89, no. B11, pp. 9305-9320, 1984.

[76] R. A. Ketcham and W. D. Carlson, "Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to the geosciences," *Computers & Geosciences*, vol. 27, no. 4, pp. 381-400, 2001.

[77] S. Iglauer, A. Paluszny, and M. J. Blunt, "Simultaneous oil recovery and residual gas storage: a pore-level analysis using in situ X-ray micro-tomography," *Fuel*, vol. 103, pp. 905-914, 2013.

[78] C. Bonnetti, M. Cuney, F. Malartre, R. Michels, X. Liu, and Y. Peng, "The Nuheting deposit, Erlian Basin, NE China: synsedimentary to diagenetic uranium mineralization," *Ore Geology Reviews*, vol. 69, pp. 118-139, 2015.

[79] T. Ohtani, Y. Nakashima, and H. Muraoka, "Three-dimensional miarolitic cavity distribution in the Kakkonda granite from borehole WD-1a using X-ray computerized tomography," *Engineering Geology*, vol. 56, no. 1-2, pp. 1-9, 2000.

[80] M. J. Heap, "The influence of sample geometry on the permeability of a porous sandstone. Geoscientific instrumentation," *Methods and Data Systems*, vol. 8, no. 1, pp. 55-61, 2019.

[81] I. J. Farquharson, M. J. Heap, Y. Lavallée, N. R. Varley, and P. Baud, "Evidence for the development of permeability anisotropy in lava domes and volcanic conduits," *Journal of Volcanology and Geothermal Research*, vol. 323, pp. 163-185, 2016.

[82] M. J. Heap and B. M. Kennedy, "Exploring the scale-dependent permeability of fractured andesite," *Earth and Planetary Science Letters*, vol. 447, pp. 139-150, 2016.

[83] M. J. Heap, A. R. Kushnir, H. A. Gilg, F. B. Wadsworth, T. Reuschlé, and P. Baud, "Microstructural and petrophysical properties of the Permo-Triassic sandstones (Buntsandstein) from the Soultz-sous-Forêts geothermal site (France)," *Geothermal Energy*, vol. 5, no. 1, p. 26, 2017.

[84] A. R. Kushnir, M. J. Heap, and P. Baud, "Assessing the role of fractures on the permeability of the Permo-Triassic sandstones at the Soultz-sous-Forêts (France) geothermal site," *Geothermics*, vol. 74, pp. 181-189, 2018.

[85] S. Zhou, X. Zhuang, and T. Rabczuk, "A phase-field modeling approach of fracture propagation in poroelastic media," *Engineering Geology*, vol. 240, pp. 189-203, 2018.

[86] M. O. Saar and M. Manga, "Depth dependence of permeability in the Oregon Cascades inferred from hydrogeologic, thermal, seismic, and magmatic modeling constraints," *Journal of Geophysical Research: Solid Earth*, vol. 109, no. B4, 2004.

[87] R. A. Freeze and J. A. Cherry, *Groundwater*, vol. 604, Prentice-Hall, Inc Englewood cliffs, New Jersey, 1979.

[88] L. Guillou-Frottier, C. Carré, C. Bourgine, V. Bouchot, and A. Genter, "Structure of hydrothermal convection in the Upper Rhine Graben as inferred from corrected temperature data and basin-scale numerical models," *Journal of Volcanology and Geothermal Research*, vol. 256, pp. 29-49, 2013.

[89] P. J. Linstrom and W. G. Mallard, "The NIST Chemistry WebBook: a chemical data resource on the internet," *Journal of Chemical & Engineering Data*, vol. 46, pp. 1059-1063, 2001.

[90] G. Launay, *Hydrodynamique des systèmes minéralisés périgranitiques: étude du gisement à W-Sn-(Cu) de Panasqueira (Portugal)*, PhD thesis, University of Orléans, France, 2018.

[91] L. Smith and D. S. Chapman, "On the thermal effects of groundwater flow: 1. Regional scale systems," *Journal of Geophysical Research: Solid Earth*, vol. 88, no. B1, pp. 593-608, 1983.

[92] M. Rabinowicz, J. Boulegue, and P. Genthon, "Two- and three-dimensional modeling of hydrothermal convection in the sedimented Middle Valley segment, Juan de Fuca ridge," *Journal of Geophysical Research: Solid Earth*, vol. 103, no. B10, pp. 24045-24065, 1998.

[93] C. Garibaldi, L. Guillou-Frottier, J. M. Lardeaux et al., "Thermal anomalies and geological structures in the Provence basin: implications for hydrothermal circulations at depth," *Bulletin de la Société Géologique de France*, vol. 181, no. 4, pp. 363-376, 2010.

[94] C. Martel, R. Champallier, G. Prouteau et al., "Trachyte phase relations and implication for magma storage conditions in the Chaine des Puys (French Massif Central)," *Journal of Petrology*, vol. 54, no. 6, pp. 1071-1107, 2013.

[95] F. H. Cornet and D. Burlet, "Stress field determinations in France by hydraulic tests in boreholes," *Journal of Geophysical Research: Solid Earth*, vol. 97, no. B8, pp. 11829-11849, 1992.

[96] Y. Gunzburger and V. Magnenet, "Stress inversion and basement-cover stress transmission across weak layers in the Paris basin, France," *Tectonophysics*, vol. 617, pp. 44-57, 2014.

[97] K. Reiter and O. Heidbach, "3-D geomechanical-numerical model of the contemporary crustal stress state in the Alberta
Basin (Canada)," *Solid Earth*, vol. 5, no. 2, pp. 1123-1149, 2014.

[98] T. J. Buchmann and P. T. Connolly, "Contemporary kinematics of the Upper Rhine Graben: a 3D finite element approach," *Global and Planetary Change*, vol. 58, no. 1-4, pp. 287-309, 2007.

[99] D. P. Yale, "Fault and stress magnitude controls on variations in the orientation of in situ stress," *Geological Society, London, Special Publications*, vol. 209, no. 1, pp. 55-64, 2003.

[100] D. L. Whitney and B. W. Evans, "Abbreviations for names of rock-forming minerals," *American Mineralogist*, vol. 95, no. 1, pp. 185-187, 2009.

[101] D. Cassard, J.-C. Chabod, E. Marcoux, B. Bourgine, and C. Castaing, "Mise en place et origine des minéralisations du gisement filonien à Zn, Ge, Ag, (Pb, Cd) de Noailhac-Saint-Salvy (Tarn, France)," *Chronique de la recherche minière*, vol. 514, pp. 3-37, 1994.

[102] M. Jébrak, *Contribution à l'histoire naturelle des filons F-Ba des Hercynides françaises et marocaines*. Document BRGM No. 99, p. 510, 1984.

[103] M. Jébrak, "Hydrothermal breccias in vein-type ore deposits: a review of mechanisms, morphology and size distribution," *Ore Geology Reviews*, vol. 12, no. 3, pp. 111-134, 1997.

[104] W. J. Phillips, "Hydraulic fracturing and mineralization," *Journal of the Geological Society*, vol. 128, no. 4, pp. 337-359, 1972.

[105] M. Jébrak, "Les textures intra-filoniennes, marqueurs des conditions hydrauliques et tectoniques," *Chronique de la Recherche Minière*, vol. 506, pp. 25-35, 1992.

[106] S. G. Hagemann, D. I. Groves, J. R. Riddley, and J. R. Vearn-combe, "The Archean lode gold deposits at Wiluna, Western Australia; high-level brittle-style mineralization in a strike-slip regime," *Economic Geology*, vol. 87, no. 4, pp. 1022-1053, 1992.

[107] Y. Nara, P. G. Meredith, T. Yoneda, and K. Katsuhiko, "Influence of macro-fractures and micro-fractures on permeability and elastic wave velocities in basalt at elevated pressure," *Tectonophysics*, vol. 503, no. 1-2, pp. 52-59, 2011.

[108] G. H. Eggertsson, Y. Lavallée, J. E. Kendrick, and S. H. Mar-kússon, "Improving fluid flow in geothermal reservoirs by thermal and mechanical stimulation: the case of Krafla vol-cano, Iceland," *Journal of Volcanology and Geothermal Research*, vol. 391, p. 106351, 2020.

[109] P. Sruoga, N. Rubinstein, and G. Hinterwimmer, "Porosity and permeability in volcanic rocks: a case study on the Serie Tobifera, South Patagonia, Argentina," *Journal of Volcanology and Geothermal Research*, vol. 132, no. 1, pp. 31-43, 2004.

[110] M. J. Heap, B. M. Kennedy, J. I. Farquharson et al., "A multi-disciplinary approach to quantify the permeability of the Whakaari/White Island volcanic hydrothermal system (Taupo Volcanic Zone, New Zealand)," *Journal of Volcanology and Geothermal Research*, vol. 332, pp. 88-108, 2017.

[111] J. L. Cant, P. A. Siratovich, J. W. Cole, M. C. Villeneuve, and B. M. Kennedy, "Matrix permeability of reservoir rocks, Ngatamariki geothermal field, Taupo Volcanic Zone, New Zealand," *Geothermal Energy*, vol. 6, no. 1, p. 2, 2018.

[112] M. J. Heap, V. R. Troll, A. R. Kushnir et al., "Hydrothermal alteration of andesitic lava domes can lead to explosive volcanic behaviour," *Nature Communications*, vol. 10, no. 1, pp. 1-10, 2019.

[113] M. J. Heap, D. M. Gravley, B. M. Kennedy, H. A. Gilg, E. Bertolett, and S. L. Barker, "Quantifying the role of hydrothermal alteration in creating geothermal and epithermal mineral resources: the Ohakuri ignimbrite (Taupō Volcanic Zone, New Zealand)," *Journal of Volcanology and Geothermal Research*, vol. 390, article 106703, 2020.

[114] J. Elder, *Geothermal systems*, Academic Press, London, 1981.

[115] A. Bejan and S. Lorente, "The constructual law and the evolution of design in nature," *Physics of Life Reviews*, vol. 8, no. 3, pp. 209-240, 2011.

[116] L. Gosselin and A. Bejan, "Constructual heat trees at micro and nanoscales," *Journal of Applied Physics*, vol. 96, no. 10, pp. 5852-5859, 2004.

[117] R. N. Horne, "Three-dimensional natural convection in a confined porous medium heated from below," *Journal of Fluid Mechanics*, vol. 92, no. 4, pp. 751-766, 1979.

[118] A. Morris, D. A. Ferrill, and D. B. Henderson, "Slip-tendency analysis and fault reactivation," *Geology*, vol. 24, no. 3, pp. 275-278, 1996.

[119] D. A. Ferrill, J. Winterle, G. Wittmeyer et al., "Stressed rock strains groundwater at Yucca Mountain, Nevada," *GSA Today*, vol. 9, no. 5, pp. 1-8, 1999.

[120] D. L. Siler and J. E. Faulds, *Three-dimensional geothermal fairway mapping: examples from the western Great Basin, USA (No. DOE-Pyramid-2842-dls-5)*, Geothermal Resources Council, Davis, CA (United States), 2013.

[121] C. E. Manning and S. E. Ingebritsen, "Permeability of the continental crust: implications of geothermal data and metamorphic systems," *Reviews of Geophysics*, vol. 37, no. 1, pp. 127-150, 1999.

[122] S. E. Ingebritsen and M. S. Appold, "The physical hydrogeology of ore deposits," *Economic Geology*, vol. 107, no. 4, pp. 559-584, 2012.

[123] S. E. Ingebritson and T. Gleeson, *Crustal Permeability: Introduction to the Special Issue*, 2014.

[124] O. Kolditz, H. Shao, W. Wang, and S. Bauer, *Thermo-Hydro-Mechanical Chemical Processes in Fractured Porous Media: Modelling and Benchmarking*, Springer, Berlin, 2016.

[125] S. F. Cox, "Coupling between deformation, fluid pressures, and fluid flow in ore-producing hydrothermal systems at depth in the crust," *Economic Geology*, vol. 100, pp. 39-75, 2005.

[126] R. J. Brown and J. Korringa, "On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid," *Geophysics*, vol. 40, no. 4, pp. 608-616, 1975.

[127] K. Runesson, D. Perić, and S. Sture, "Effect of pore fluid compressibility on localization in elastic-plastic porous solids under undrained conditions," *International Journal of Solids and Structures*, vol. 33, no. 10, pp. 1501-1518, 1996.

[128] H. B. Motra and H. H. Stutz, "Geomechanical rock properties using pressure and temperature dependence of elastic P- and S-wave velocities," *Geotechnical and Geological Engineering*, vol. 36, no. 6, pp. 3751-3766, 2018.