# Fluid flow in crustal fault zones with varying lengthwise thickness: application to the Margeride fault zone
(French Massif Central)

Emmy Penhoët, Laurent Arbaret, Laurent Guillou-Frottier, Hugo Duwiquet, Charles Gumiaux, Mathieu Bellanger

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

## To cite this version:

Emmy Penhoët, Laurent Arbaret, Laurent Guillou-Frottier, Hugo Duwiquet, Charles Gumiaux, et al.. Fluid flow in crustal fault zones with varying lengthwise thickness: application to the Margeride fault zone (French Massif Central). Geothermal Energy, 2025, 13, pp.6. 10.1186/s40517-025-00334-9. hal-04892277

HAL Id: hal-04892277

https://hal.science/hal-04892277v1

Submitted on 16 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Fluid flow in crustal fault zones with varying lengthwise thickness: application to the Margeride fault zone (French Massif Central)

Emmy Penhoët<fnref n="1" />,<fnref n="2" />* , Laurent Arbaret<fnref n="2" />, Laurent Guillou-Frottier<fnref n="1" />,<fnref n="2" />, Hugo Duwiquet<fnref n="3" />, Charles Gumiaux<fnref n="2" /> and Mathieu Bellanger<fnref n="4" />

## Abstract

Crustal fault zones, holding promise as potential geothermal reservoirs, remain largely untapped and unexplored. Located in the southern Massif Central, France, the Margeride fault zone (MFZ) varies in thickness (lateral extension perpendicular to the fault plane) from 100 m to over 2500 m. Reactivated several times under different stress regimes since the Variscan orogeny, this zone is characterized by an intense alteration and fracturing. As a result, the multiple reactivation of the fault zone has maintained permeability, leading to favourable conditions for fluid circulation. Structural measurements and geological cross sections were used to precisely constrain thickness and geometry of the fault zone. North of the MFZ, the Coren thermal spring indicates reservoir temperatures of about 200–250 °C, hinting at the possible existence of a temperature anomaly. To investigate this geothermal potential, 3D numerical models simulating fluid circulation within a fault zone were conducted. Various configurations were explored, altering fault zone thickness and permeability for two key geometries. The first geometry, which manipulated the width of the fault zone along its length, demonstrated a direct correlation between fault zone thickness and amplitude of thermal anomaly. Thinner faults (< 500 m) exhibited multiple weak positive thermal anomalies, while thicker faults (> 500 m) tended to develop a single, substantial positive thermal anomaly. In the second examined geometry, where fault zone thickness increased longitudinally, a consistent positive temperature anomaly emerged at the thickest section of the fault zone. Depending on the permeability value, an additional anomaly may develop but will migrate laterally towards the thinnest part of the fault zone. This multi-disciplinary approach, combining numerical modelling and field measurements, presents a predictive methodology applicable to geothermal exploration in analogous basement domains. In our case, it has shown that the northern end of the Margeride fault zone could represent an area that needs to be explored further to assert its high geothermal potential. Our numerical models will increase understanding of how fault width and geometry impact the geothermal potential of the Margeride fault zone and similar areas in crystalline basement.

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

<footnote>*Correspondence:
penhoetemmy@gmail.com</footnote>

<footnote n="1">¹ Georesources Division, BRGM,
45060 Orléans, France</footnote>

<footnote n="2">² ISTO, UMR 7327, Université d'Orléans, CNRS, BRGM,
45071 Orléans, France</footnote>

<footnote n="3">³ Geothermal Division, ENGIE,
Faubourg de l'Arche, 1 Place
Samuel de Champlain,
92930 Paris La Défense, France</footnote>

<footnote n="4">⁴ TLS-Geothermics, 91 Chemin de Gabardie, 31200 Toulouse, France</footnote>
**Keywords:** French Massif Central, Margeride fault zone, Thickness variation, 3D numerical modeling, Fault density, Structural measurements

**Introduction**

In high-temperature geothermal systems ($T > 150$ °C) where thermal convection dominates, deep hot fluids require highly permeable pathways to reach the surface (Manning and Ingebritsen 1999; Ledésert et al. 2009; Ingebritsen and Appold 2012; Ingebritsen and Gleeson 2014; Moeck 2014; Caspari et al. 2020, Jolie et al. 2021). These conduits can be faults or fractures networks or alteration pathways (Lowell 1975; Berkowitz 2002; Simms and Garven 2004, Bischoff et al. 2024) in which thermally driven convective circulation may develop. In the last decades, 2D and 3D numerical models of fluid flow in permeable fault zones have been performed to better constrain the key parameters that control the circulation of fluids and heat within the fractured system but also in the surrounding rocks (López & Smith 1995; Bächler et al. 2003; Patterson et al. 2018; Duwiquet et al. 2022). At the scale of the fault zone, most significant parameters are the permeability and thickness of the damaged zone, its horizontal and vertical extension, and inclination and orientation of the fault with respect to the regional stress field (Duwiquet et al. 2021). Indeed, many studies have shown that the faults affected by an extensional regime are excellent upwelling conduits for deep fluids (Barton et al. 1995; Mckenna and Blackwell 2004; Faulds et al. 2010; Haines et al. 2016, Roche et al. 2018). Thickness of the fault zone plays also a role in the triggering of thermal convection, as demonstrated theoretically by Malkovsky and Magri (2016). In terms of geometry, crustal fault zones that are sufficiently thick (fault zone thickness greater than 100 m) and long (length greater than 1 km) with a displacement of 100–150 m are more likely to host fluid circulation (convection), as proposed by Guillou-Frottier et al. (2024). Ultimately, topography is also a key parameter to consider in geothermal exploration, since lateral pressure gradients may control how the meteoric fluids travel through the first kilometers of the crust (Forster and Smith 1989; Barton et al. 1995; Taillefer et al. 2018, Alt-Epping et al. 2021, Carbajal-Martínez et al. 2024).

In purely crustal context, because of the time evolution of the direction and intensity of the stress field, faults zones often show multiple successive brittle reactivations. It results in the development of relatively thick fault zones whose permeability could be maintained over time. Conversely, recent faults zones that experienced a small number of brittle events without a significant change in the stress conditions may remain relatively narrow. Therefore, geothermal exploration in highly fractured crustal setting should include the identification of thick polyphased faults and the knowledge of their brittle evolution (Mitjanas et al. 2024).

To investigate the impact of the brittle evolution and consequent fault zone thickness on the geothermal potential of crustal fault zones, we selected the Margeride Fault in the French Massif Central. The French Massif Central is characterized by a large-scale surface heat flow anomaly (100–110 mW m⁻², over several hundreds of km Lucazeau and Vasseur 1989). When the radioactive component of surface heat flow is removed, the remaining mantle heat flow would be around 60 mW m⁻² (Lucazeau et al. 1984) beneath the Cenozoic grabens and the volcanic areas, implying temperatures above 200 °C at a depth of 5 km (see Fig. 5 in Lucazeau et al. 1984). Numerous crustal fault zones or
permeable fault clusters exist in the Massif Central (Pontgibaud fault zone, Duwiquet et al. 2019; border faults of the tertiary Limagne Basin, Calcagno et al. 2014; Chaudes-Aigues fault network, Vasseur et al. 1997). In the Saint Flour area (Fig. 1), the Margeride Fault is the result of a long-lived brittle history since the Carboniferous. Consequently,

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Geological setting of the Margeride fault zone (MFZ). Coordinate system: Lambert 93—RGF93 EPSG:2154. a location of the study area. b Simplified geological map of the Margeride Granitic Complex issued from the 1:1 000 000 geological map of mainland France (BRGM, https://www.brgm.fr/en/news/news/open-science-free-access-brgm-digital-geological-maps). c Simplified map of the study area. Blue lines locate the four selected sections broadly oriented perpendicular to the MFZ and crossing the four major half grabens. Blue dots refer to outcrops on which statistical structural measurements of fractures orientations where performed and that are located as close as possible to the geological sections. Green dots refer to representative locations (Loc, see text for details). d Close-up map of the northwestern segment of the MFZ along which the Cenozoic Saint Flour Basin emplaced</figcaption></figure>
the Margeride fault is made of a complex network of gouge alternating with damage zones that we refer hereafter as the Margeride fault zone (MFZ). Furthermore, the MFZ has the particularity of having an increasing thickness from south-east to north-west.

Aside from the structural geology arguments, a thermo-mineral spring is located in the thickest part of the MFZ. The "Font de Vie" spring at Coren-les-Eaux, (10 km north of Saint-Flour) is found in the gneissic bedrock. Even if the spring is cold at the surface (10 °C), it is nevertheless highly mineralized (7 g L⁻¹), and geothermometers indicate at depth a reservoir temperature of 220 °C (Fouillac 1983).

We performed structural measurements of the regional brittle structures network in the granitic and gneissic surrounding rocks and in the Margeride fault zone itself. We particularly focused along four sections striking perpendicularly to the MFZ and centered on the associated cenozoic grabens. Results have been used to constrain the 3D numerical models, where the variation in thickness from south-east to north-west is taken into account.

Given the few studies on the effects of varied fault widths on fluid and heat transfer (see, however, Duwiquet et al. 2024), our work based on key points, such as the brittle evolution of the fault zone, and on numerical models will help to constrain the identification of new potential geothermal targets.

## Geological setting

The regional-scale crustal fault zone used as a study-case consists of the Margeride fault zone (named MFZ hereafter) located in the central part of the French Massif Central (Fig. 1). The French Massif Central is essentially inherited from the Variscan collisional orogen whose successive ductile to brittle orogenic events range from the late Silurian (ca. 410 Ma) to the early Permian (ca. 300 Ma, Faure 1995; Vanderhaeghe et al. 2020). In addition, from the late Paleozoic to present-day, this granite and gneiss basement has suffered successive tectonic brittle events (Blès et al. 1989; Faure et al. 1995; Merle et al. 2023).

The MFZ crosscuts several variscan lithologies: a late variscan granite (named the Margeride Granite Complex hereafter) and its country rock composed of two Paleozoic gneissic series. In addition, Oligocene half grabens developed at the contact of the MFZ evidence its brittle cenozoic reactivation (Fig. 1).

## The metamorphic units

In the Saint Flour-Le Malzieu area (Fig. 1c), the metamorphic formations related to the Devonian to early Carboniferous, south verging, ductile shearing are represented by the Lower and Upper Gneiss units (named, respectively, LGU and UGU hereafter; Vanderhaeghe et al. 2020; Fig. 1). The UGU is made of retrogressed granulite composed by garnet-bearing diatexites associated with metatexites. The LGU consists mainly in micaschistes, paragneisses and orthogneisses that suffered greenschist to amphibolitic metamorphism. The contact between the two units is interpreted as a main shear zone between two stacked nappes (Burg and Matte 1978).
## The Margeride granitic complex (MGC)

The 323 ± 12 Ma old Margeride granite complex (MGC, red area in Fig. 1c) emplaced during a late carboniferous extensional event as a laccolith within the LGU and UGU series (Couturié 1979; Talbot et al. 2005). The main facies are composed of a porphyritic monzogranite crosscut by late leucogranites and quartz dykes. The emplacement model proposed by Talbot et al. (2005) involves four NE–SW trending feeding zones that played either as pull-apart or extensional gashes controlled by sinistral NW–SE trending faults. Mineral lineations developed during the magmatic emplacement are homogenously oriented SE–NW with a shallow dip to the NW (Talbot et al. 2005) it implying a syn-orogenic emplacement within an NW–SE extensional regional stress field (Vanderhaeghe et al. 2020).

The MGC shows several signs favorable to the circulation of fluids. There are numerous quartz, porphyrites and microgranites dykes which could be potential pathways for fluids. Evidence come from the reactivation of N140E brittle fractures that allow the upwelling of basaltic melts during the Cenozoic [Julianges plateau described by Couturié (1972), Fig. 1b]. Fifteen km to the West of Saint Flour (Fig. 1c), the granite has been suggested as a hydrothermal reservoir associated with the Chaudes-Aigues thermal springs, based on the above criteria (Cailleaux et al. 1976; Vasseur et al. 1997).

## The Margeride fault zone (MFZ)

The MFZ is a major crustal fault zone whose age and multiphase history remain poorly documented. Broadly, the MFZ could be divided into two segments: The northern segment that begins south of the Saint Flour Basin (Fig. 1) is composed of a ~N130E trending, 15-km-long and up to 2500-m-thick network of brittle faults that overprint a similarly oriented, preexisting ductile fault zone affecting the LGU and UGU series. Its northwestward extension is occulted by the basaltic lava flows of the mio-pliocene Cantal stratovolcano. This structure is easily visible in the topography, since there is an important slope break between the Margeride plateau culminating at an elevation of 1400 m to the East and the Malzieu Basin at 800 m to the West. The age of the early ductile event is not determined. However, its orientation, parallel to the mineral lineations in the UGU and LGU series, together with the fact that it does not propagate to the south into the Margeride granite point out that this ductile event ended before 320 Ma. It could, therefore, be attributed to the first major ductile extension of the late Variscan period (Faure 1995).

The second southeastern segment transects the Margeride granite. It is composed of a ~N150E trending brittle fault system. The thickness of the MFZ is not observable in the field. However, the minimum distance between the outcropping continental sediments and the unfractured granite visible along the Malzieu and Saint-Alban Basins highlights a thickness of 100 m, at most. The lengthwise connection with the northern segment takes place near the Chaliers Basin, where the transition from N160E to N140E trending along both faults and grabens is obvious (see Fig. 1, Burg and De Goer de Herve 1990).
## Late Variscan and Mesozoic brittle overprint

During the late Variscan period (330–280 Ma), the region experienced successive major extensions. First, an NW–SE extension related to the syn-collisional extension that affected the ductile part of the crust. Second, an NE–SW extension (Faure 1995; Daffos et al. 2024). The latest variscan-related event was a North–South compression that progressively became East–West during the Stephanian (Blés 1989). This North–South compression generated a conjugate brittle system made of N140E–N170E dextral and N010E–N050E sinistral strike-slip faults, associated with N80–N120 inverse faults. At the end of the Paleozoic, the Variscan range was peneplaned.

Tectonic events associated with the opening of the Atlantic and Ligurian oceans are understood to have little effect on the study area, as evidenced by an accentuated peneplain process during the Mesozoic.

## Cenozoic extension and uplift

In the late Cretaceous and Paleocene, an intraplate compression involved a buckling of the lithosphere and reactivation of pre-existing crustal discontinuities (Ziegler et al. 2007). During the Eocene, an east–west extensional regime developed over the Massif Central with the establishment of the European Cenozoic Rift System (ECRIS; Michon and Merle 2001; Dèzes et al. 2004; Merle et al. 2023). Although the origin of the ECRIS remains debated (Michon and Merle 2001; Michon et al. 2003; 2005; Dèzes et al. 2004; 2005), all authors agree that a mantle thermal anomaly started during the Paleocene and caused thinning of the lithosphere beneath the FMC. It resulted in the implementation of N140E to N0E elongated grabens, from North to South of the FMC, the northernmost being the Limagne Basins.

In the study area, four half grabens are located along the western side of the MFZ. From the South to the North, they comprise the Saint-Alban Basin, the Malzieu Basin, the Chaliers Basin and the Saint-Flour Basin (De Goer De Herve et al. 1990; Fig. 1). Sedimentation was predominantly detrital. All Basins are limited on their eastern edge by the MFZ, along which subsidence was limited as the thickness of the continental sediments has been estimated to a maximum of 150 m (De Goer De Herve et al. 1990). The western edges of these half-grabens are not marked by major faults, leading to their asymmetric geometry. In the Malzieu Basin, this asymmetric structure is outlined by the present course of the Truyère River (not shown in Fig. 1).

During the Upper Miocene, the Massif Central was subjected to a regional uplift due to an NW–SE directed compression (Mazabraud et al. 2005). It resulted in the active erosion of the peneplaned surface leading to the formation of the Truyère valley.

## Present-day regional stress field

According to Mazabraud et al. (2005) the current stress field in the French Central Massif implies an NW–SE compression, which means that fractures oriented in this direction could play as tensions gashes. Indeed, favourably oriented faults will be those that enable fluid circulation (Barton et al. 1995; Mckenna and Blackwell 2004,
Faulds et al. 2010, Duwiquet et al. 2022). Thus, the two segments of the MFZ, which varies in strike by 30°on average, will possibly react differently in this regional stress field.

## Geothermal setting of the Malzieu area

### Surface and mantle heat flow

Along the MFZ, several surface heat flow values have been published (Lucazeau and Vasseur 1989). They are consistent as a heat flow value of 96 mW m⁻² is found in the North (Saint Flour area, site “Saint Michel”, NGHF Database, Lucazeau 2019), while a heat flow of 102 mW m⁻² is obtained in the south (“Aumont Aubrac” site, Lucazeau 2019). A geothermal gradient of 31 °C km⁻¹ has been measured for the first site and 34 °C km⁻¹ for the second. However, these sites correspond to shallow boreholes (depth of 150 m), where measured temperatures could be affected by paleoclimatic disturbances. Consequently, we chose two others neighboring heat flow sites in the FMC (Auriat and Villechabrolle), with depths of 1000 and 650 m, respectively. These two boreholes show deep unperturbed temperature gradients of 30 and 32 °C km⁻¹, respectively. These values are consistent with the average temperature gradient of the FMC.

The large-scale surface heat flow anomaly of the FMC has been interpreted by Lucazeau et al. (1984) as a mantle diapir bringing an additional mantle heat flow of 25–30 mW m⁻², thus leading to a total mantle heat flow around 60 mW m⁻² in the area.

### Radiogenic heat production

Decades ago, four boreholes were drilled by the INAG (National Institute for Astronomy and Geophysics) in the MGC to collect core samples and to measure the temperature gradient. Heat production measurements on these core samples indicate values varying between 2.76 µW m⁻³ and 4.32 µW m⁻³. The value finally retained is 3.75 µW m⁻³ (Bayer et al. 1982). For the metamorphic units, samples taken over the FMC point to a mean value of 2.5 µW m⁻³ (Lucazeau and Vasseur 1981).

## Methods

### Brittle structures analysis

To constrain the brittle structures affecting the MGC and its surrounding gneissic basement outcrop-scale measurements at key sites were combined with regional scale fault density analysis methods, as detailed below.

### Outcrop-scale structural measurements

We performed outcrop scale fractures analysis on 57 references sites, 31 of them being distributed along, and as close as possible to, the four selected cross sections (Fig. 1c). These sections are perpendicular to, and crosscut, the MFZ. In addition, they are centered on the four cenozoic half-grabens. In each location, from 50 to 200 measurements were performed according to the size of the outcrops. Low dipping joints, i.e., dipping less than 15°, were not considered in these locations as they more likely relate to either exfoliation or unloading joints. Nevertheless, two additional stations were selected to characterize the change in attitude of these flat-lying joints approaching Le Malzieu graben.
Due to the absence of adequately sized and representative outcrops, measurements were not conducted within the continental sedimentary series filling the grabens.

## Regional scale fault density

A regional scale faults density analysis has been performed on an assemblage composed of the four geological maps covering the study area, namely, N°789-Saint-Flour (De Goër de Herve et al. 1990), N°790-Langeac (Bouiller et al. 1986), N°813-Chaudes-Aigues (Burg et al. 1991), N°814-Saugues (Couturié 1979). We digitized faults by applying a grid consisting of 500 m wide square cells. The cumulative length values per cell were divided into 5 classes using the natural breaks of the distribution histogram (Jenks algorithm implemented in QGIS; Jenks 1967).

## Primary mineral structures analysis

In addition to the analysis of brittle structure, mineral foliation and lineation developed in the two gneissic units were measured at 57 stations distributed over the area. Foliations and lineations in the MGC are issued from Talbot et al. (2005). Talbot et al. (2005) evidenced the close correspondence between the magnetic fabrics and the mineral fabrics in the granite. Both magmatic and ductile foliations were used to properly define the fractures type developed in the MGC and the surrounding basement, respectively, as they may overlay on primary joints (Marre 1986). In igneous rocks, three main type of primary joints, which are orthogonal to each other and guided by the orientation of inherited mineral foliation and lineation, may develop (Cloos 1922): the Q joints (cross fractures) that are perpendicular to the mineral lineation, the S-joints (longitudinal fractures) that are perpendicular to the foliation and parallel to the lineation, and finally, the L joints that are parallel to the foliation (Marre 1986).

## Numerical modeling

The thickening of the MFZ, from the Le Malzieu to the Saint-Flour areas may have an impact on the possibility for thermal convection to occur within the permeable fault zone. Before investigating by numerical modeling whether thermal convection can occur in various parts of the Margeride fault zone, a simple analytical approach is presented below.

### Estimate of Rayleigh number

Thermal convection occurs if the Rayleigh number $Ra$ (ratio of driving versus resisting forces, depending on physical properties, medium geometry, and temperature contrast) exceeds a critical value $Ra_c$ (e.g., Turcotte and Schubert 2002, chap. 9.9). For an infinite porous layer with constant temperature conditions on top and bottom boundaries, this $Ra_c$ value equals $4\pi^2$ (Lapwood 1948).

For a fault zone of a finite width, the critical value ($Ra_c$) has been studied by Malkovsky and Magri (2016) through a linear stability analysis and for a temperature-dependent fluid viscosity. For a fault zone with a half-width over height ratio $\Delta$, they showed that for an averaged temperature gradient, the critical Rayleigh number writes:
$$Ra_c = 2.466 \left[ \left( \frac{8.19}{\Delta} \right)^{5/4} + (4\pi^2)^{5/4} \right]^{4/5} \quad (1)$$

In our case, the average temperature gradient hypothesis is valid, since values around 30 °C km⁻¹ have been measured (see above). For a constant viscosity case, the coefficient 2.466 in Eq. (1) is replaced by 1 (see Malkovsky and Magri 2016), and the classic 4π² value is recovered as soon as the half-width over height ratio (Δ) is large (i.e., for an infinite width). For a finite fault zone thickness, the critical Rayleigh number is larger than 4π².

Assuming that the Margeride fault zone height is 5 km (which corresponds to the gneiss-granite transition according to Talbot et al. 2005), fault zone half-width can be estimated at 50 m for the Malzieu area, while it reaches 1250 m in the Saint-Flour area. The associated Δ values are thus 0.01 and 0.25, respectively. It follows that the two critical Rayleigh numbers for Le Malzieu and Saint-Flour areas equal:

$$Ra_c(\text{Malzieu}) = 2056$$

$$Ra_c(\text{Saint-Flour}) = 155$$

In other words, the thick deformation zone in the Saint-Flour area may easily host thermal convection, since the critical Rayleigh number is small. For given physical conditions and petrophysical properties, one may observe thermal convection in the Saint Flour deformation zone, while these conditions will not lead to thermal convection in the Le Malzieu area.

In the following, different fault zone geometry (thicknesses and lengths) and different physical properties will be tested through numerical modeling of heat and mass transfer within a finite-width permeable fault zone. Numerical models were performed with the Comsol Multiphysics™ software.

### 3D models

As demonstrated by Malkovsky and Magri (2016), the establishment of a convective regime depends directly on the morphology of the fault zone and more precisely on its thickness and height. A parametric study where fault zone thickness and fault zone permeability are varied is first presented. To study the impact of the varying thickness in the longitudinal direction, fault zone with a prismatic shape is then studied.

As in Guillou-Frottier et al. (2013) Darcy's law was coupled with heat transfer equation and mass conservation. Several variables such as thermal conductivity, heat capacity,

<table><caption>**Table 1** Parameters used for rocks and fluid in 3D modelling</caption><thead><tr><th></th><th>Porosity</th><th>Permeability K<sub>0</sub> (m<sup>2</sup>)</th><th>Matrix thermal conductivity (W/m.K)</th><th>Density (kg/m<sup>3</sup>)</th><th>Matrix heat capacity (J/kg.K)</th><th>Heat production (W/m<sup>3</sup>)</th><th>Fluid thermal conductivity (W/m.K)</th><th>Fluid heat capacity (J/kg.K)</th></tr></thead><tbody><tr><td>Fault zone</td><td>0.20</td><td>10<sup>-15</sup> - 5 10<sup>-14</sup></td><td>3.30</td><td rowspan="2">2700</td><td rowspan="2">1000</td><td>x</td><td rowspan="2">0.6</td><td rowspan="2">4200</td></tr><tr><td>Bulk</td><td>0.05</td><td>10<sup>-18</sup></td><td>2.80</td><td>3 10<sup>-6</sup></td></tr></tbody></table>

These values have been obtained from different studies (Smith and Chapman 1983, Rabinowicz et al. 1998, McKenna and Blackwell 2004, Guillou-Frottier et al. 2013, Taillefer 2017, Duwiquet et al. 2019)
porosity, rock density (see Table 1) necessary for modeling were estimated thanks to field observations and literature data.

**Boundary and initial conditions** At the upper surface, a temperature of 10 °C was fixed. According to regional data, a constant heat flux of 85 mW m⁻² is imposed at the base of the model (Lucazeau and Vasseur 1989). An averaged heat production rate of 3 µW m⁻³ (see previous section on radiogenic heat production rate) is assumed in the entire domain The lateral borders were thermally insulated (host rocks). A geothermal gradient of 30 °C km⁻¹ was imposed as an initial condition. The equivalent thermal conductivity is deduced from the above thermal conditions (heat flow, temperature gradient and heat production), leading to a value of 2.8 W m⁻¹ K⁻¹, consistent with thermal conductivity measurements in the area (Lucazeau et al. 1984).

The surface pressure is equivalent to the atmospheric pressure (10⁵ Pa), and its value increases gradually with depth. A no flow condition was imposed on lateral sides and bottom of the model (host rocks).

**Meshing and geometry** For the parametric study, a simplified parallelepidic geometry is chosen for the fault zone and the basement (Fig. 2a). The size of mesh elements is based on the permeability value. In areas of low permeability (host rocks), the mesh element size is coarse (minimum size 100 m) but extremely fine in the fault zone (high permeability) (maximum 20 m). The same mesh features have been used for the prismatic shape.

**Variable physical properties** For the fluid, here supposed as pure water, the dynamic viscosity (Pa.s) is a function of temperature T(°C), such that:

$$ \mu(T) = 2.414 \cdot 10^{-5} \cdot e^{\frac{570}{T+133}} \qquad (2) $$

following Smith and Chapman (1983) and Rabinowicz et al. (1998).

The density of water is also a function of pressure and temperature, such that:

<figure><img src="image_6.png" /><figcaption>**Fig. 2** a Simplified geometry used for parametric studies, with coarse mesh in the embeddings and refined mesh in the faulted zone; the fault zone thickness h varies between 100 and 2500 m. b Models with a lengthwise variation of the fault zone thickness from 2500 to 50 m (prismatic shape). The length of the fault zone, L, has been varied, but the meshes shown here correspond to the results of Figs. 12 and 13 where L = 20 km, close to the MFZ. The height of the fault zone is 4 km, and it crops out the surface</figcaption></figure>
$$ 
\begin{aligned}
\rho_f(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 (3) $$

Equation (3) is based on a polynomial fit of experimental data, as developed by Launay (2018) and published in Guillou-Frottier et al. (2020).

The main unknown in the modelling approach remains the permeability value for each formation found in our study area. The MGC and the two metamorphic units are considered, in the absence of fracturing, as impermeable. The permeability values in common rocks vary between $10^{-18}$ m² and $10^{-13}$ m² (Brace 1984; Evans et al. 1997; Sonney and Vuataz 2009; Ranjram et al. 2015; Walter 2016; Achtziger-Zupancic et al. 2016; 2017; Neuzil 2019), which leaves a wide range of possibilities. Here we will attempt to better constrain qualitatively the permeability using structural measurements, field observations and fault density estimates. The increase in pressure at depth tends to seal the fractures and thus causes a decrease in the permeability of the units. This phenomenon can be described by the following equation:

$$ K_u(z) = K_{0u} e^{-\frac{z}{\delta}} \quad (4) $$

where $K_u(z)$ is the permeability of each unit (m²), $K_{0u}$ is the highest permeability (from $10^{-15}$ m² to $5 \times 10^{-14}$ m² in the fault zone, and $10^{-18}$ m² in the host rocks) reached by the unit (at the surface) and z is the depth (m) (Garibaldi et al. 2010). The length δ (m) characterizes the intensity of the decrease in permeability with depth. A δ value of 2000m has been chosen.

The following numerical results correspond to final state (steady state or pseudo steady state regimes) of 2 Myr of transient evolution. Initial conditions reproduce a conductive regime and at t > 0, permeabilities are effective.
# Results

## Structures analyses

### Primary ductile and magmatic structures

In the two gneissic units the mineral foliations show a similar distribution marked by two plane populations. The main population corresponds to an N065 trending mean plane that dips to the North-West by 33° (Fig. 3a). The secondary population yields a mean foliation oriented N105E and dipping to the North by 64°. In contrast the lineation measured on these two associated foliation populations yield a mean lineation constantly trending N100E and plunging to the Northwest by 20° (Fig. 3b). This N100E regional stretching lineation is known to be developed over the southern metamorphic nappes of the French Massif Central (Faure 1995). These two associated foliations, characterized by different dips but close strike directions represent flanks of asymmetric folds developed in both gneissic units.

These ductile structures of Eo-Variscan age are characterized by a first N100 direction, which is the regional stretching lineation and is also widely developed in the southern metamorphic nappes of the Massif Central (Faure 1995).

<figure><img src="image_7.png" /><figcaption>**Fig. 3** Equal area projection (lower hemisphere) of primary ductile and magmatic mineral fabrics. **a** Foliations and **b** lineations measured in 57 locations distributed in both UGU and LGU units. Kamb contours in standard deviation. **c** AMS foliations and **d** lineations of the MGC (modified from Talbot et al. 2005). All stereographic projections are in an equal area lower hemisphere</figcaption></figure>
In the MGC, field measurements of mineral foliation and lineation are difficult to achieve, as the granite appears macroscopically isotropic (no visible magmatic foliation). Nevertheless, Talbot et al. (2005) measured the magnetic foliations and lineations by applying the anisotropy of magnetic susceptibility technic and they considered them as closely representative of the magmatic mineral foliation and lineation, respectively. Talbot et al. (2005) evidenced magnetic foliations weakly dipping either to the NE or the SW, leading to a ~N45E trending pole girdle (Fig. 3c). In contrast, the tied magnetic lineations are well organized as they are almost horizontal with an average trending direction of N130E to N160E (Fig. 3d).

## Regional fracture analysis

Regional fracture network analysis was performed by combining 57 locations distributed in the study area. For a sake of simplification and as most of the fractures exhibit a strong dip, strike directions are also presented in rose diagrams together with lower hemisphere projections (Fig. 4a, b). Among these locations, three reference sites have been selected as they best illustrate the fracture network in both UGU and LGU gneissic units, the MGC and finally, in a location close to the MFC (Fig. 4c; locations in Fig. 1).

In the LGU and UGU units, 941 fractures were measured (Fig. 4a). The main population is represented by N20E to N30E trending, almost vertical, fractures. A second population is composed of well oriented N90E trending planes. Two minor, more distributed, N45E and N140E trending populations complete this fracture network. The main population is well represented in the selected reference site Loc1 (Fig. 4c), in which the second population trend N70E.

In the MGC, 3345 fractures were measured (Fig. 4b). The main population is composed of vertical fractures trending from N140E to N160E. The second population is even more scattered as it is formed of vertical fractures with strikes ranging from N45E

<figure><img src="image_8.png" /><figcaption>**Fig. 4** a Azimuthal (left) and orientation (right) distributions of fractures measured in the UGU and LGU gneissic units. b Azimuthal (left) and orientation (right) distribution of fractures measured in the MGC. c orientation distributions in three reference sites in the UGU (Loc1), in the MGC (Loc2) and along the MFZ (Loc3, see Fig. 1 for geographical locations). All stereographic projections are in an equal area lower hemisphere. All Kamb contours in standard deviation</figcaption></figure>
to N80E. This fracture network organization is well expressed in the reference site Loc2 (Fig. 4c).

Approaching the MFZ, the fracture network changes (Loc3, Fig. 4c). Only one population represented by well-organized N150E to N160E trending fractures is found. This fracture network is parallel to the MFC at this location (Fig. 1).

## Flat lying joints

Flat lying, joints are commonly observed in the MGC. These joints could be related to either L joints or unloading joints. A reference site outcrop West of Saint-Alban-Sur-Limagnole (Loc 4, Fig. 5a; see location in Fig. 1c). In this outcrop, a main vertical, N120E trending, fracture population is crosscut by, regularly distributed, flat-lying joints (Fig. 5b). Approaching the grabens, and more particularly near Le Malzieu graben, the attitude of the joints changes. This is particularly notable in the Porte des Fées Canyon, where the attitude of the joints is clearly visible in the landscape (Loc 5, Fig. 5c; see location in Fig. 1c). In this location two vertical, N170E and N60E trending fracture populations are crosscut by N160E trending joints gently dipping to the East of 12° (Loc 5, Fig. 5d).

<figure><img src="image_9.png" /><figcaption>**Fig. 5** a roadside outcrop located West of Saint-Alban-Sur-Limagnole (Loc 4, GPS coordinates: 44.7787, 3.3537, see location on Fig. 1d) showing the flat-lying joints in the MGC. b Orientation distribution of primary fractures and joints in Loc 4. c Outcrop located on the western edge of the Malzieu basin (Loc 5, GPS coordinates: 44.8632, 3.3171, see location on Fig. 1d) showing tilted joints in the MGC. d Orientation distribution of primary fractures and joints in Loc 5</figcaption></figure>
Cross sections

Four cross sections perpendicular to the MFZ and centered over the four sedimentary Basins were designed. Two are located in the MGC to the South (Saint-Alban sur-Limagneole and Chinchazes Basins; Le Malzieu Basin; Fig. 6) and two are located in the LGC/UGC gnessic basement to the north (Chaliers Basin; Saint Flour Basin; Fig. 7) (locations in Fig. 1c). No outcrop appropriated for statistical fracture analysis were found in the oligocene sediments of the grabens.

<figure><img src="image_10.png" /><figcaption>(b)</figcaption></figure>

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

**Fig. 6** a Geological section crossing the Saint-Alban-sur-Limagnole Basin and Chinchazes Basin. Loc 4 refers to flat-lying joints analysis location (see Fig. 5a, b). b Geological section crossing Le Malzieu Basin. For (a) and (b), both scales are in meters. Loc 5 refers to flat lying joints analysis location (see Fig. 5a, b). All stereographic projections are in an equal area lower hemisphere. All Kamb contours in standard deviation. Principal stress tensors were calculated using FaultKin software
<figure><img src="image_12.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 7** a Geological section crossing the Chaliers Basin. b Geological section crossing the Saint Flour Basin. All stereographic projections are in an equal area lower hemisphere. For (a) and (b), both scales are in meters. All Kamb contours in standard deviation. Principal stress tensor was calculated using FaultKin software</figcaption></figure>

Seven reference sites are distributed near the section crossing over the two Basins of Saint-Alban-sur-Limagnole and Chinchazes (Fig. 6a). Among them, site 41 was dedicated to flat lying joints analysis (Loc 4, Fig. 5a, b). Two main populations of fractures are recognized. Northeast of the Basins, an N050E-to-N070E trending family is associated with an N100E family. Both are vertical in average. Stress tensor analysis in location 49 indicates that a horizontal East–West main regional stress field was responsible for this conjugate fault system (Fig. 6a). These fracture families are present westward and around the sedimentary Basin, while they progressively
rotate to form N110E and N170E trending, strongly dipping, fractures population in the westernmost location (location 50, Fig. 6a). Secondary fracture family comprise N170E to North–South planes that were only observed in locations 40 and 43.

Eight references sites document the fracture distribution along the section cross-
ing over Le Malzieu Basin (Fig. 6b). Among them, site 37 was dedicated to flat-lying
joint analysis (Loc 5, Fig. 5c, d). The two main vertical fractures family trending
N050E to N070E and N100E to N150E, respectively, are present all along the sec-
tion. Stress tensor analysis in sites 63 and 52 indicate that a horizontal NW–SE main
regional stress field was responsible for this conjugate fault system (Fig. 6a). Sec-
ondary fractures comprise N020E to North–South vertical planes (sites 51 and 55),
while N090E trending vertical fractures are observed in sites 54 and 52 (Fig. 6b).

The MFZ is not outcropping along the two sections. Therefore, the thickness can
only be estimated using indirect evidence. Along the eastern border of the Chin-
chazes Basin, the undeformed granite outcrop less than hundred meters from the
oligocene continental sediments filling the graben (outcrop location: 44.7826,
3.4158). This apparent low thickness is consistent to those proposed for the segment
of the MFZ bordering Le Malzieu Basin (Couturié 1972).

Further to the North, the Chaliers Basin is emplaced in the LGU-UGU units (Fig. 7a). 7 reference sites are distributed along the section. To the East of the MFZ bordering this graben the gneisses are affected by North-South to N030E highly dipping fractures only associated occasionally with a secondary, almost East-West trending, fracture population (Sites 18 and 17; Fig. 7a). West of the graben, near Chaliers village, these two fracture populations are also present. Conjugate N010E and N160E fault system was found near the contact with the continental sediments filling the graben (Site 24, Fig. 7a). Stress tensor analysis yields a horizontal NW-SE main regional stress field similar to the one obtained near Le Malzieu graben (Fig. 6b), all being compatible with a broadly N070E-directed horizontal extension.

The northern last section crosses over the Saint Flour Basin (Fig. 7b). East of
the graben, the LGU outcrops along the road of Montchamp (Margeride plateau,
Fig. 8a). From the East to the West, the gneiss that exhibits an N120E trending, 80°
West dipping, foliation (Fig. 8b), is progressively affected by an increasing ductile
deformation that becomes more pervasive with evidence of mylonitisation while
approaching the border of the graben. This ductile deformation is crosscut by dis-
tributed brittle normal faults almost parallel to the ductile foliation (Fig. 8c). 8 refer-
ence sites for fracture analysis are distributed along the section. To the East, fracture
analysis evidence two populations formed of N010E to N030E fractures associated
with N100E fractures, both being vertical (Sites SF7, Fig. 7b). These two popula-
tions are also present westward, where they progressively rotate from 20° clockwise
approaching the border of the graben. This zone of about 500 m thick is composed
of isolated lenses of intensively fractured gneiss lenses and quartz veins embedded
in an argilised breccia (Fig. 9). Within the graben, these two populations are pre-
sent but from place to place, one take precedence over the other. In location SF4
(Fig. 7b), N040E trending fracture is the main population, while to the West, the
N120E to N140E trending fractures become the only population (location SF5 and
SF8, Fig. 7b).
<figure><img src="image_14.png" /><figcaption>**Fig. 8** a Roadside section of the eastern border of the MFZ, East of the Saint Flour graben, depicting the transition from the fresh gneiss with a strongly westward dipping foliation (b, 50m away from the view) to a damaged and altered zones (Loc6, GPS coordinates: 45.0534, 3.1797, location in Fig. 1d). c Youngest normal, west dipping, faulting event crosscutting the locally mylonitized gneiss</figcaption></figure>

## Regional fault density

The fault density distribution complements the outcrop scale fracture analyses performed along the four main geological sections by evidencing regional scale major fault zones.

In the MGC, low to moderate fault densities domains form N140E trending faults closely parallel to the MFZ and distributed over the granitic massif (Fig. 10). An exception comes from the Chaudes-Aigues area, 20 km. South-West of our study area, where the fault density is very high forming also an N140E trending lineament but that propagate northwestward in the LGU near Chaudes-Aigues, a city known for its hot springs swarm (De Goër de Herve and Burg 1991). It should be noticed that the fault density remains low along the trace of the MFZ, which is an additional argument for a thin core zone forming this segment emplaced in the MGC.

The second domain comprising moderate to high fault density is formed by the N045E trending contact between the MGC and the LGU (section between Chaudes-Aigues area and Le Malzieu Basin, Fig. 10). On the contrary, East of the MFZ, the contact between the MGC and the basement is not marked by a significant fault density.

In the North-East, the gneissic units show several N045E trending lineaments of moderately fault densities.

Finally, the larger domain formed of very high fault density of faults corresponds to the N120E directed segment of the MFZ bordering the Saint Flour Basin (Fig. 10). Its maximum width is about 2500 m.
<figure><img src="image_15.png" /><figcaption>**Fig. 9** a Global and b zoomed-in field views of the fault core of the MFZ, East of the Saint Flour graben, composed of isolated lenses of fractured gneiss lenses embedded in an altered breccia (Loc7, GPS coordinates: 45.0533, 3.1604, location in Fig. 1d)</figcaption></figure>

## Synthesis

By combining field observation and fracture network analysis at outcrop scale with regional fault density analysis we can draw up the structural constraints of the
<figure><img src="image_16.png" /><figcaption>**Fig. 10** Fault density map of the study area, the contact between lithologies (black and red lines) correspond to those shown in Fig. 1. Coordinate system: Lambert 93—RGF93 EPSG:2154</figcaption></figure>

magmatic and tectonic history of the region since the Variscan period that led to the lengthwise change in thickness and orientation of the MFZ.

The southeastern segment of the MFZ is emplaced in the MGC and is composed of a, N140Etrending, narrow fault zone of possibly less than 100m thick (this study, Couturié 1972). This orientation, found over the magmatic body in the form of distributed lineaments of high fault density (Fig. 10), corresponds also to the direction of the four main feeding zones that opened as transtensional sinistral relays and that lead to the emplacement of the 310–300 Ma old MGC (Talbot et al. 2005). Magmatic fabrics recorded these directions of emplacement with N150E trending lineations on roughly flat foliations (Fig. 3c, d). Consequently, longitudinal primary joints are vertical and strike N140E to N150E. Such longitudinal joints were reactivated, shortly after the cooling of the main granitic facies and until 290 Ma, leading to the successive emplacement of microgranites, then microporphyrite, and finally, quartz dykes (De Goër de Herve and Burg 1991). Stress tensor analysis on conjugate fault system is compatible with a horizontal NW–SE principal stretching axis (Fig. 6). Such direction prevailed during the late carboniferous contemporaneously with dykes’ emplacement, but also during the Pliocene, a period
during which the region experienced an intense uplift and a subsequent erosion (Ziegler and Dèzes 2007; Daffos et al. 2024).

The northwestern segment of the MFZ is emplaced in the UGU and LGU and is composed of a, N120E trending, up to 2500 thick fault zone (this study, De Goër de Herve and Tempier 1988). De Goër de Herve and Tempier (1988) showed that an N140E trending and vertical, muscovite-rich, mylonitic shear zone affected the gneisses of the two units, before the Visean, i.e., before the emplacement of the MGC. This vertical, sinistral strike-slip fault, continued to operate as a brittle fault, particularly during the Stephanian; leading to a thick fault zone formed of densely fractured and argilised preexisting mylonites (Fig. 9, De Goër de Herve and Tempier 1988).

All these variscan brittle structures were reactivated during the two successive active extension events during the Oligocene and then the late Miocene (De Goer De Herve et al. 1990). Reactivation as normal faults of both N140E trending primary joints and similarly oriented late variscan faults in the granite together with reactivation of the thick fault zone in the LGU-UGU led to the subsidence of the four half-grabens. These normal movements along west dipping faults, to the east, and evidence of quasi stratigraphic unconformable contact on the underlying granite, along western edges. To explain this geometry, Couturié (1979) and De Goer De Herve et al. 1990 proposed a listric geometry of the eastern boundary faults of the Basins, which are themselves subsiding in a "Roll-Over" fashion (Figs. 6 and 7). Change in attitude of the flat lying joints near le Malzieu graben is an additional evidence for a "Roll-Over" rotation of about 12° (Fig. 5). Finally, during the Mio-Pliocene compressive recession (Merle and Michon 2001), N40 trending fractures played as vertical relays leading to the present-day "piano key" geometry of the grabens along the MFZ (Couturié (1979), De Goer De Herve et al. 1990).

The high density of faults (Fig. 10) and the various reactivations of fragile structures indicate that the MFZ could have a high permeability value. Several studies have shown that the permeability of fault zones can reach values of the order of $10^{-13} \text{ m}^2$ and ever more (Neuzil et al. 2019, Scibek 2020). In our study we will, therefore, use a range of values from $10^{-15} \text{ m}^2$ to $5 \ 10^{-14} \text{ m}^2$.

The MFZ, that comprises a southern N140E trending, thin fault zone, branching to the northwest to a 2500 m thick fault zone represents a study case of crustal fault zone with lengthwise change in thickness. The influence of this geometry on the thermal regime and upward fluid transfer is approached by numerical modeling.

## Numerical modeling

The entire numerical procedure has been previously tested with several benchmark experiments, including the 3D benchmark tests of Magri et al. (2017), as detailed in Guillou-Frottier et al. (2020).

## Sensitivity analysis

A series of numerical experiments has been performed with different fault zone thicknesses (100 to 2500 m) and different fault zone permeabilities. Temperature field, positive thermal anomalies (defined as the excess temperature when compared to the conductive regime), velocity field, and convective patterns, are recorded.
The value of the maximal thermal anomaly depends on the thickness of the fault zone (Fig. 11a). The optimum thickness (1000 m) for achieving the highest thermal anomaly value (55 °C) is not the largest, since the value reached for a fault zone thickness of 2500 m is 52 °C. However, if we compare the temperature anomalies produced at the same depth (1000m), this time it is the thickest fault zones that generate the highest values: 48 °C for 2500 m compared with 35 °C for 500 m (see dashed red line in Fig. 11a). The impact of thickness is directly observed. For a fixed permeability value (K), for example at $K=3 \ 10^{-14} \ m^2$ (star symbol in Fig. 11a), the temperature anomaly values produced are not the same. Thickest faults, over 1000 m, produce higher anomaly values than thinner faults. At low permeabilities ($7 \ 10^{-15} \ m^2$, circle symbol in Fig. 11a), the anomaly increases with thickness. On the contrary, for higher permeabilities ($5 \ 10^{-14} \ m^2$), an optimum thickness of 1000 m appears. For thicker or thinner faults, the anomalies decrease. Thick fault zones will be more favorable for deep anomalies generating higher temperatures than superficial anomalies in thinner fault zones.

The depth of the 150 °C isotherm also varies with the fault zone thickness (Fig. 11b). For low permeability ($5 \ 10^{-15} \ m^2$), heat conduction dominates, and the depth is common to all models. When the permeability value increases, two distinct behaviors emerge for thin faults (100–1000 m) and thick faults (1700–2500 m). In the first case, the depth of the 150 °C isotherm decreases as permeability increases, reaching 3800 m for a fault zone thickness of 500 m. In the second case, convection operates, cold downwellings take over upwellings and causes this boundary to deepen below 4000 m. It appears that a thick permeable zone favors the downward flow of cold fluids and causes the fault zone to cool. The variation in thickness shows that for fine faults (<1000 m), a rise in the 150°C isotherm will lead to a high increase in the permeability of the fault zone ($3 \ 10^{-14} \ m^2$). The same phenomenon can be observed for thick faults, but with permeability values an order of magnitude lower.

<figure><img src="image_17.png" /><figcaption>**Fig. 11** Graphs showing the impact of variations in the permeability K and thickness of the fault zone on the value of the maximum thermal anomaly (a), the depth of the 150 °C isotherm (b) and fluid flow velocity (c)</figcaption></figure>
For a better understanding of the different thermal regimes, the fluid flow velocities have been recorded (Fig. 11c). Velocity values depend on the width of the fault zone. The greater the thickness, the higher the maximum velocity. When $h=2500$ m, velocities reach $9 \times 10^{-8}$ m s$^{-1}$ compared to $10^{-11}$ m s$^{-1}$ for 100 m or $2 \times 10^{-8}$ m s$^{-1}$ for 500 and 1000 m. Two behaviors again stand out in terms of the increase in these velocities as permeability increases. The finer fault zones (500–100 m) show a curve that can be divided into two areas: a first part with a small increase from $10^{-15}$ m$^2$ to $7 \times 10^{-15}$ m$^2$ and a second part with a higher slope up to $5 \times 10^{-14}$ m$^2$. For thicknesses greater than 500 m, the curve splits into several portions and leads to higher values. Once again, reasoning with a fixed permeability directly highlights the role of thickness variation. For a permeability of $10^{-14}$ m$^2$, a 500 m fault zone has a velocity of $10^{-11}$ m/s compared with $5 \times 10^{-8}$ m/s for a 2500 m fault zone. This difference, of three orders of magnitude, indicates that at a fixed permeability, the thickness of the fault zone is crucial for the establishment of a convective regime.

The results of the various models can be divided into one conductive and two convective regimes. The first group represents the configurations leading to a conductive regime (green area, in Fig. 12). The associated fluid flow velocities are very low, around $10^{-12}$ m s$^{-1}$, and isotherms remain horizontal (no temperature distortion).

When permeability increases, fluid velocities also increases, and thermal convection is triggered. The first convective regime corresponds to a steady-state single upwelling (II, yellow area, in Fig. 12). As described in other studies (e.g., Zhao et al. 2003; Guillou-Frottier et al. 2020) the upwelling zones exhibit a triangular-like pattern. Fluids rise in the center of the structure and are bordered by two large downwellings on each side.

<figure><img src="image_18.png" /><figcaption>**Fig. 12** a Three convective regimes as a function of fault zone thickness and permeability, black dots represent all the models carried out to constrain the regime boundaries (dashed black lines). b Three examples of the different regimes as defined in (a). All these experiments have been performed with a fault zone length L equal to 20 km</figcaption></figure>
Fluid circulation velocities are higher than in the previous case, between 2.5 $10^{-9}$ m s$^{-1}$ (500 m thickness) and $9 \times 10^{-8}$ m s$^{-1}$ (2500 m thickness). Positive maximum thermal anomalies range from 38 °C (500 m) to 55 °C (1000 m).

The second convective regime presents several convection cells evidenced by several zones of upwelling (III, orange area, in Fig. 12). This regime has already been highlighted in other studies. As permeability increases, the wavelength between the upwelling zones decreases (Guillou-Frottier et al. 2020). In our case, two upwellings appear and move laterally, along the longitudinal direction (transient regime). In the case illustrated on bottom-right of Fig. 12 (fault zone thickness of 500 m and permeability of $3 \times 10^{-14}$ m$^2$), two upwelling zones are visible, bordered by downward flows. Here again, the upwelling exhibits a triangular like pattern, as described in Guillou-Frottier et al. (2020). Since permeability decreases with depth, cold dense fluids move more easily from the most permeable zone (surface) and become predominant over hot upwellings, which start to rise in the less permeable zone (bottom). This phenomenon gives them a triangular shape, with upwellings edges “eroded” by cold downwellings. One consequence is that the generated temperature anomalies have lower values. Fluid flow velocities vary between $10^{-8}$ m s$^{-1}$ and $9 \times 10^{-8}$ m s$^{-1}$.

The boundaries between the different domains are not horizontal lines, indicating that these regimes are directly dependent on the thickness of the fault zone. The boundary between the conductive regime and the first convective regime decreases with thickness (between blue and yellow areas, Fig. 12a). At a fixed permeability, the development of convection is favored when the thickness of the permeable zone increases.

For fine fault zone (<500 m), the establishment of a type III convective regime appears to be directly affected by the variation in thickness. For a permeability of $5 \times 10^{-14}$ m$^2$, a fault 100 m thick will be in a conductive regime, one 250 m thick in convective regime II and one 500 m thick in regime III. However, beyond 500 m, the limit seems to be guided more by permeability. As the thickness increases, the permeability value that allows the third regime to be reached also increases. For a fixed permeability value, the limit becomes almost horizontal (between orange and yellow areas, Fig. 12a). A wider range of stability is observed, for thickness greater than 500 m, leading to single steady-state upwelling (yellow area, see Fig. 12a).

For fault zones less than 250 m thick, conduction is the dominant regime (for realistic permeabilities). Indeed, for thin fault zones, convection occurs only for high permeability (e.g., for $h=100$ m, $K=5 \times 10^{-13}$ m$^2$) in transient regime, so that the third regime is favored. Furthermore, a stationary system is difficult to achieve, which leads directly to the third regime (orange area, Fig. 12a).

## Models with a lengthwise variation of the fault zone thickness

Field evidence as well as fault density distribution (Fig. 10) underline the north-westward thickening of the MFZ. Numerical models of Fig. 13 include a lengthwise varying thickness, from 70 m (bottom right of the fault zone) to 2500m (upper left). The maximum fault zone permeability $K_{F0}$ equals $8 \times 10^{-15}$ m$^2$ in Fig. 13a, and $1.5 \times 10^{-14}$ m$^2$ in Fig. 13b. In addition to temperature field illustrated as in Fig. 12, positive thermal anomalies $\Delta T$ are represented by isosurfaces separated by 10 °C (see colour bars). In both experiments, there is no thermal anomalies in the thinnest parts of the fault zone, where
<figure><img src="image_19.png" /><figcaption>**Fig. 13** Results of the models with a lengthwise variation of the fault zone thickness for two permeabilities</figcaption></figure>

isotherms flatten (conductive regime). On the opposite, a thermal anomaly is always present around the thickest part of the fault zone (convective regime). Its amplitude can vary with time but not its location. Experiment of Fig. 13a reaches a steady-state 1.5 Myr after the beginning of the experiment. On the opposite, experiment of Fig. 13b exhibits a quasi-periodic behaviour, where thermal anomalies migrate laterally towards the thinnest zones (black arrows). At $t_0 + 500$ kyr, three thermal anomalies are visible, including the fixed one around the thickest part of the fault zone. With time, the central one decreases in amplitude and disappears, while the third one increase in amplitude, as shown at time $t_0 + 2$ Myr. The dashed black arrows indicate the new incipient upwelling, which will also migrate toward the thinnest zone. Top views on the right of Fig. 13b illustrate the thermal diffusion of the anomaly within the embeddings.

# Discussion

## Exploration using a multidisciplinary approach

The aim of this work was to study the effects of the progressive thickening of a simple fault towards a thick fault zone on the circulation of fluids at depth. Structural measurements in the field were used to constrain qualitatively the intensity of fracturing and the thickness of the fault zone. These data have been used directly in the 3D dynamic models to demonstrate the direct impact of thickness on convective dynamics. Our numerical results indicate that deep anomaly values (1000 m depth) increase with thickness (Fig. 11c). Furthermore, thick fault zones (> 500 m) are more likely to host stable convective systems. Finally, the latest models have demonstrated the existence of a positive temperature anomaly (Fig. 13), whatever the permeability used, at the thickest end of the fault zone. To achieve full coupling between field exploration and dynamic modeling, the measured orientations would still have to be used in a model incorporating the
current stress regime of the study area. Such models, which are beyond the scope of this study, need a thermo-hydro-mechanical coupling approach (e.g., Duwiquet et al. 2021, 2022).

## Controlling fluid flow factors depend on the thickness of the fault zone

As expected, the higher permeability leads to an increase in the number of upwelling in the fault zone (Rabinowicz 1998, Guillou-Frottier et al. 2020). The first boundary shown in Fig. 12 between the conduction domain and a single upwelling (blue and yellow areas) seems logical, since the increase in width coupled with the increase in permeability is favourable to the development of convection (Magri Malkovsky 2016). Presented as the parameter governing fluid flow, permeability does not appear to be the only controlling parameter for the regimes shown in Fig. 12. In fact, as the width of the fault zone increases, so does the yellow area represent stable convection. Indeed, it appears that a thick fault zone allows the development of steady-state convection. We might, therefore, assume that fault zones less than 500 m thick are mainly governed by permeability, but that this is no longer true for thicker fault zones.

## Edge effects and transient regime

Concerning the transient aspect of the third regime (zone in red, Fig. 12), this migration of upwelling zones appears to be due to the dimensions of the fault zone. Indeed, different lengths of the fault zones have been tested ($L = 15, 20, 50$ km) and edge effects have been observed. When fault zone length is compatible with the convective wavelength associated with a given permeability (e.g., Guillou-Frottier et al. 2020, their Fig. 1), then a steady-state regime establishes. If the length of the permeable zone and the wavelength of the upwelling zones are not compatible, the convective pattern is not stable and the upwelling migrates laterally towards the thinnest parts of the fault zone, then disappearing and followed by a new incipient upwelling around the center of the fault zone, which will migrate with time (dashed arrow in Fig. 13b). It is, however, important to point out that whatever the convective regime, the thermal anomaly located at the thickest part of the fault zone does not migrate and remains fixed in space and time.

## Geothermal potential of the Margeride fault zone

The advantage coupling field observations and modelling (define the brittle of the study area to determine the state of permeability and ultimately use thermal modelling) is that it can be used to provide assistance before the exploration phase, making it possible to predict areas with high geothermal potential. The use of numerical models has already confirmed the possibility to identify exploitable areas independently of other exploration methods, as for the Rittershoffen area, close to the Soultz-sous-Forêts geothermal site (Upper Rhine graben, France) (Guillou-Frottier et al. 2013). In our case we have shown that thick fault zones are more likely to develop stable convection leading to high temperature anomalies (Fig. 12). Models with a lengthwise thickness variation indicate that a positive temperature anomaly persisted at the thick end of the fault zone. The thin end is not affected by convection movements. Following these results, a positive temperature anomaly of 50 to 70 °C could be expected in the thickest part of the fault zone, or in the intermediate part at a depth of around 1000 m. These models tend to show that
high-temperature geothermal exploitation would be unsuitable in view of the temperatures reached (between 100 and 120 °C at a depth of 1000 m). However, the results obtained here are encouraging, given that a temperature of 10 °C is imposed at the surface of the model (see Boundary Conditions section), which could considerably minimise the value of the anomaly due to the imposed cooling.

Studies such as those by Eldursi et al. (2023) and Duwiquet et al. (2022) have also shown that the orientation of the fault zone fundamentally influences fluid circulation. According to Mazabraud (2005), the zone is currently undergoing NW-SE compression, leading to the opening of similarly oriented structures. The northern part of the MFZ (Saint-Flour), therefore, appears to be favorably oriented (N120E on average), unlike the southern part (Malzieu, N150E). Finally, the presence of the Font-de-Vie Spring at Coren-les-Eaux is an important indicator, since it indicates a reservoir at a depth of over 220 °C (Fouillac 1983).

Our numerical approach applied to the MFZ should be improved by more realistic (stress regime, boundary condition, topography) modelling. Nevertheless, it indicates that geothermal exploration should focus primarily on the north of the study area (Saint-Flour), while the south of the zone shows a very low potential (Malzieu).

## Conclusions

The MFZ is a natural laboratory, enabling us to study a range of processes that have a direct impact on fluid circulation. Using this example, it was possible to quantify the impact of the thickness of a fault zone on fluid circulation and convective regimes. Thin fault zones will be more likely to produce multiple upwellings with weaker, shallower anomaly values. At the opposite, it will be more common to obtain a single higher and deeper anomaly in thick fault zones. It, therefore, appears that in addition to permeability, the thickness of the fault zone becomes a fundamental parameter to be taken into account. Finally, we showed that, for fault zones with longitudinal thickening, the thickest extremity will always be affected by a positive temperature anomaly.

This method has identified a potential area for more extensive exploration in crystal-line basement. However, some elements still need to be included in dynamic models. The current stress field is not taken into account, nor is the change in orientation of the MFZ. Finally, the fault zone is marked by a major jump in topography (around 600 m) that could influence fluid circulation. These new elements will help to refine the geothermal potential of the study area.

### Abbreviations

MFZ Margeride fault zone

BRGM Bureau de Recherches Géologiques et Minières

ISTO Institut des Sciences de la Terre d'Orléans

FMC French Massif Central

LGU Lower Gneiss Unit

UGU Upper Gneiss Unit

MGC Margeride Granite Complex

ECRIS European Cenozoic Rift System

### Author contributions

This work corresponds to the first year of PhD thesis of EP. LA collected, analysed and interpreted the field data. LGF worked on numerical models. They both actively participated in writing the manuscript. HD helped analysed the results of numerical models. CG collected, analysed field data. MB took part in a field mission to collect data. All authors read and approved the final manuscript.
**Funding**

The financial support was provided by Bureau de Recherches Géologiques et Minières (BRGM), Institut des Sciences de la Terre d'Orléans (ISTO) and by the project « GERESFAULT » (ANR-19-CE05-0043).

**Availability of data and materials**

All data are available from the corresponding author on reasonable request.

**Declarations**

**Competing interests**

The authors declare that they have no competing interests.

Received: 27 April 2024 Accepted: 18 December 2024

Published online: 16 January 2025

**References**

Achtziger-Zupančić P, Loew S, Hiller A, Mariethoz G. 3D fluid flow in fault zones of crystalline basement rocks (Poehl-Tellerhaeuser Ore field, Ore Mountains, Germany). Geofluids. 2016;16(4):688-710.

Achtziger-Zupančić P, Loew S, Mariéthoz G. A new global database to improve predictions of permeability distribution in crystalline rocks at site scale. J Geophys Res Solid Earth. 2017;122(5):3513-39.

Alt-Epping P, Diamond LW, Wanner C, Hammond GE. Effect of glacial/interglacial recharge conditions on flow of meteoric water through deep orogenic faults: insights into the geothermal system at Grimsel Pass, Switzerland. J Geophys Res Solid Earth. 2021;126(7):e2020JB021271.

Bächler DL. Impact of graben-parallel faults on hydrothermal convection—Rhine Graben case study. Phys Chem Earth Parts A/B/C. 2003;28(911):431-41.

Barton CA, Zoback MD, Moos D. Fluid flow along potentially active faults in crystalline rock. Geology. 1995;23:683-668.

Bayer R, Couturie JP, Vasseur G. Résultats géophysiques récents sur le granite de la Margeride. Annales De Géophys. 1982;38:431-47.

Berkowitz B. Characterizing flow and transport in fractured geological media: a review. Adv Water Resour. 2002;25(8-12):861-84.

Bischoff A, Heap MJ, Mikkola P, Kuva J, Reuschlé T, Jolis EM, et al. Hydrothermally altered shear zones: a new reservoir play for the expansion of deep geothermal exploration in crystalline settings. Geothermics. 2024;1(118): 102895.

Blés JL, Bonijoly D, Castaing C, Gros Y. Successive post-Variscan stress fields in the French Massif Central and its borders (Western European plate): comparison with geodynamic data. Tectonophysics. 1989;169(1):79-111.

Bouiller R, Cornen G, Burg JP, Lasnier B. Carte géologique de la France à 1/50 000, feuille 790 Langeac. BRGM; 1986.

Brace WF. Permeability of crystalline rocks: New in situ measurements. J Geophys Res. 1984;89(B6):4327-30.

Burg JP, De Goer de Herve A. Carte géologique de la France à 1/50 000, feuille 813 Chaudes-Aigues. BRGM; 1991.

Burg J-P, Matte P. A cross section through the French Massif central and the scope of its Variscan geodynamic evolution. Z Dtsch Geol Ges. 1978;129:429-60.

Cailleaux P, Fouillac C, Michard G, Ouzounian G. Etude géochimiques des sources thermales de Chaudes-Aigues (Cantal). Conséquences Goéthermiques C R Acad Sci. 1976;282(13):1237-40.

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

Carbajal-Martínez D, Wanner C, Diamond LW, Peiffer L, Fletcher JM, Inguaggiato C, et al. Behavior of amagmatic orogenic geothermal systems: insights from the Agua Blanca fault, Baja California, Mexico. Geochem Geophys Geosyst. 2024;25(3):e2023GC011145.

Caspari E, Greenwood A, Baron L, Egli D, Toschini E, Hu K, et al. Characteristics of a fracture network surrounding a hydrothermally altered shear zone from geophysical borehole logs. Solid Earth. 2020;11(3):829-54.

Cloos E. "Feather joints" as indicators of the direction of movements on faults, thrusts, joints and magmatic contacts. Geology. 1932;18(5):387-95.

Couturié JP. Carte géologique détaillée de la France à 1/50 000, feuille 814 Saugues. BRGM; 1972.

Couturié J-P, Vachette Caen M, Vialette Y. Age Namurien d'un laccolite granitique différencié par gravité : Le granite de la Margeride (Massif Central Français). Acad Sci. 1979;289(5):449-52.

Daffos C, Arbaret L, Bourdier J-L, Gumiaux C. Regional structural control on the Mont-Dore Plio-quaternary volcanism (France). J Volcanol Geotherm Res. 2024;447: 108029.

De Goer de Herve A, Tempier P, Boudon G. Carte géologique de la France à 1/50 000, feuille 789 Saint-Flour. BRGM; 1990.

Dèzes P, Schmid SM, Ziegler PA. Evolution of the European Cenozoic Rift System: interaction of the Alpine and Pyrenean orogens with their foreland lithosphere. Tectonophysics. 2004;389(1-2):1-33.

Duwiquet H, Arbaret L, Guillou-Frottier L, Heap MJ, Bellanger M. On the geothermal potential of crustal fault zones: a case study from the Pontgibaud area (French Massif Central, France). Geotherm Energy. 2019;7(1):33.

Duwiquet H, Guillou-Frottier L, Arbaret L, Bellanger M, Guillon T, Heap MJ. Crustal fault zones (CFZ) as geothermal power systems: a preliminary 3D THM model constrained by a multidisciplinary approach. Geofluids. 2021;2021:1-24.

Duwiquet H, Magri F, Lopez S, Guillon T, Arbaret L, Bellanger M, et al. Tectonic regime as a control factor for crustal fault zone (CFZ) geothermal reservoir in an amagmatic system: a 3D dynamic numerical modeling approach. Nat Res Res. 2022;31(6):3155-72.
Duwiquet H, Genter A, Guillou-Frottier L, Donzé FV, Ledru P, Magri F, et al. Advanced 3D TH and THM modeling to shed light on thermal convection in fault zones with varying thicknesses. *J Geophys Res Solid Earth*. 2024;129(4):e2023JB028205.

Eldursi K, Scholtes L, Conin M, Golfier F, Ledru P, Mercadier J, et al. 3-D hydromechanical simulation of intersecting faults: influences on fluid circulation and formation of oriented-orebodies. *J Struct Geol*. 2023;171: 104864.

Evans JP, Forster CB, Goddard JV. Permeability of fault-related rocks, and implications for hydraulic structure of fault zones. *J Struct Geol*. 1997;19(11):1393-404.

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

Faure M. Late orogenic carboniferous extensions in the Variscan French Massif Central, Tectonics. 1995;14(1):132-53.

Forster C, Smith L. The influence of groundwater flow on thermal regimes in mountainous terrain: a model study. *J Geophys Res.* 1989;94(B7):9439-51.

Fouillac C. Chemical geothermometry in CO2-rich thermal waters. Example of the French massif central. Geothermics. 1983;12(2-3):149-60.

Garibaldi C, Guillou-Frottier L, Lardeaux J-M, Bonté D, Lopez S, Bouchot V, et al. Thermal anomalies and geological struc- tures in the Provence basin: Implications for hydrothermal circulations at depth. *Bulletin De La Société Géologique De France*. 2010;181(4):363-76.

Guillou-Frottier L, Carré C, Bourgine B, Bouchot V, Genter A. Structure of hydrothermal convection in the Upper Rhine Graben as inferred from corrected temperature data and basin-scale numerical models. *J Volcanol Geoth Res*. 2013;256:29-49.

Guillou-Frottier L, Duwiquet H, Launay G, Taillefer A, Roche V, Link G. On the morphology and amplitude of 2D and 3D thermal anomalies induced by buoyancy-driven flow within and around fault zones. Solid Earth. 2020;11(4):1571-95.

Guillou-Frottier L, Milesi G, Roche V, Duwiquet H, Taillefer A. Heat flow, thermal anomalies, tectonic regimes and high- temperature geothermal systems in fault zones. *C R Geosc*. 2024;356(S2):1-33.

Haines S, Lynch E, Mulch A, Valley JW, Van Der Pluijm B. Meteoric fluid infiltration in crustal-scale normal fault systems as indicated by δ¹⁸O and δ²H geochemistry and ⁴⁰Ar/³⁹Ar dating of neoformed clays in brittle fault rocks. Litho-sphere. 2016;8(6):587-600.

Ingebritsen SE, Appold MS. The physical hydrogeology of ore deposits. *Econ Geol*. 2012;107(4):559-84.

Jenks GF. The data model concept in statistical mapping. *Int Yearbook Cartogr*. 1967;7:186-90.

Jolie E, Scott S, Faulds J, Chambefort I, Axelsson G, Gutiérrez-Negrín LC, et al. Geological controls on geothermal resources for power generation. *Nat Rev Earth Environ*. 2021;2(5):324-39.

Lapwood ER. Convection of a fluid in a porous medium. *Math Proc Camb Phil Soc*. 1948;44(4):508-21.

Launay G. Hydrodynamique des systèmes minéralisés péri-granitiques: étude du gisement à W-Sn-(Cu) de Panasqueira (Portugal). Université d'Orléans; 2018.

Ledésert B, Hebert R, Genter A, Bartier D, Clauer N, Grall C. Fractures, hydrothermal alterations and permeability in the Soultz Enhanced Geothermal system. *Comptes Rendus Géoscience*. 2010;342(7-8):607-15.

López DL, Smith L. Fluid flow in fault zones: analysis of the interplay of convective circulation and topographically driven groundwater flow. *Water Resour Res*. 1995;31(6):1489-503.

Lowell RP. Circulation in fractures, hot springs, and convective heat transport on mid-ocean ridge crests. *Geophys J Int*. 1975;40(3):351-65.

Lucazeau F. Analysis and mapping of an updated terrestrial heat flow data set. *Geochem Geophys Geosyst*. 2019;20(8):4001-24.

Lucazeau F, Vasseur G. Production de chaleur et régime thermique de la croûte du Massif Central. *Production de chaleur et régime thermique de la croûte du Massif Central*. 1981;

Lucazeau F, Vasseur G. Heat flow density data from France and surrounding margins. *Tectonophysics*. 1989;164(2):251-8.

Lucazeau F, Vasseur G, Bayer R. Interpretation of heat flow data in the French Massif Central. *Tectonophysics*. 1984;103(1):99-119.

Malkovsky VI, Magri F. Thermal convection of temperature-dependent viscous fluids within three-dimensional faulted geothermal systems: estimation from linear and numerical analyses. *Water Resour Res*. 2016;52(4):2855-67.

Manning CE, Ingebritsen SE. Permeability of the continental crust: implications of geothermal data and metamorphic systems. *Rev Geophys*. 1999;37(1):127-50.

Marre J. The structural analysis of granitic rocks. London: North Oxford Academic; 1986.

Mazabraud Y, Béthoux N, Guilbert J, Bellier O. Evidence for short-scale stress field variations within intraplate central-western France: intraplate short-scale stress field variations. *Geophys J Int*. 2004;160(1):161-78.

McKenna JR, Blackwell DD. Numerical modeling of transient basin and range extensional geothermal systems. *Geothermics*. 2004;33(4):457-76.

Merle O, Aumar C, Labazuy P, Merciecca C, Buvat S. Structuration tertiaire et quaternaire du Plateau des Dômes (Chaîne des Puys, Massif central, France). Géologie de la France. 2023;(1).

Merle O, Michon L. The formation of the West European Rift; a new model as exemplified by the Massif Central area. *Bull Soc Geol Fr*. 2001;172(2):213-21.

Michon L, Merle O. The evolution of the Massif Central rift: spatio-temporal distribution of the volcanism. *Bull Soc Geol Fr*. 2001;172:201-11.

Mitjanas G, Walsh JJ, Roca E, Álías G, Queralt P, Ledo J, et al. The importance of structural complexity in the localization of geothermal systems: a case study along the Vallès-Penedès fault in the Catalan coastal ranges (NE Spain). *Geothermics*. 2024;116: 102855.

Moeck IS. Catalog of geothermal play types based on geologic controls. *Renew Sustain Energy Rev*. 2014;37:867-82.

Neuzil CE. Permeability of Clays and Shales. Annu Rev Earth Planet Sci. 2019;47(1):247-73.

Patterson JW, Driesner T, Matthai SK. Self-organizing fluid convection patterns in an en Echelon fault array. *Geophys Res Lett*. 2018;45(10):4799-808.
Rabinowicz M, Boulègue J, Genthon P. Two- and three-dimensional modeling of hydrothermal convection in the sedi-
mented middle valley segment, Juan de Fuca Ridge. J Geophys Res. 1998;103(B10):24045–65.

Ranjram M, Gleeson T, Luijendijk E. Is the permeability of crystalline rock in the shallow crust related to depth, lithology or tectonic setting? Geofluids. 2015;15(1–2):106–19.

Scibek J. Multidisciplinary database of permeability of fault zones and surrounding protolith rocks at world-wide sites.
London: Nature Publishing Group; 2020.

Simms MA, Garven G. Thermal convection in faulted extensional sedimentary basins: theoretical results from finite-
element modeling. Geofluids. 2004;4(2):109–30.

Smith L, Chapman DS. On the thermal effects of groundwater flow: 1. Regional scale systems. J Geophys Res.
1983;88(B1):593–608.

Sonney R, Vuataz F-D. Numerical modelling of Alpine deep flow systems: a management and prediction tool for an exploited geothermal reservoir (Lavey-les-Bains, Switzerland). Hydrogeol J. 2009;17(3):601–16.

Taillefer A, Guillou-Frottier L, Soliva R, Magri F, Lopez S, Courrioux G, et al. Topographic and faults control of hydrothermal circulation along dormant faults in an Orogen. Geochem Geophys Geosyst. 2018;19(12):4972–95.

Talbot J-Y, Faure M, Chen Y, Martelet G. Pull-apart emplacement of the Margeride granitic complex (French Massif Central) implications for the late evolution of the Variscan orogen. J Struct Geol. 2005;27(9):1610–29.

Turcotte DL, Schubert G. Geodynamics. Cambridge: Cambridge University Press; 2002.

Vanderhaeghe O, Laurent O, Gardien V, Moyen J-F, Gébelin A, Chelle-Michou C, et al. Flow of partially molten crust con-
trolling construction, growth and collapse of the Variscan orogenic belt: the geologic record of the French Massif
Central. BSGF Earth Sci Bull. 2020;191:25.

Vasseur G, Michard G, Fouillac C. Contraintes sur la structure profonde et le fonctionnement du système hydrothermal de Chaudes-Aigues (France). 1997;(4):3–17.

Walter B. Réservoirs de socle en contexte extensif: genèse, géométries et circulations de fluides: exemples du rift intrac-
ontinental du lac Albert (Ouganda) et de la marge proximale d'Ifni (Maroc). Université de Lorraine; 2017.

Zhao C, Hobbs BE, Mühlhaus HB, Ord A, Lin G. Convective instability of 3-D fluid-saturated geological fault zones heated from below. Geophys J Int. 2003;155(1):213–20.

Ziegler PA, Dèzes P. Cenozoic uplift of Variscan Massifs in the Alpine foreland: timing and controlling mechanisms. Global Planet Change. 2007;58(1–4):237–69.

**Publisher’s Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.