# Tectonic Regime as a Control Factor for Crustal Fault Zone (CFZ) Geothermal Reservoir in an Amagmatic System: A 3D Dynamic Numerical Modeling Approach

H. Duwiquet, F. Magri, S. Lopez, T. Guillon, L. Arbaret, M. Bellanger, L.

Guillou-Frottier

► To cite this version:

H. Duwiquet, F. Magri, S. Lopez, T. Guillon, L. Arbaret, 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. Natural Resources Research, 2022, 31 (6), pp.3155-3172. 10.1007/s11053-022-10116-w. insu-03836300

HAL Id: insu-03836300

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

Submitted on 4 Nov 2022

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
<figure><img src="image_2.png" /></figure>

# Tectonic Regime as a Control Factor for Crustal Fault Zone (CFZ) Geothermal Reservoir in an Amagmatic System: A 3D Dynamic Numerical Modeling Approach

H. Duwiquet<sup>1,2,3,6,7</sup> F. Magri<sup>4,5</sup> S. Lopez,<sup>1</sup> T. Guillon,<sup>1</sup> L. Arbaret<sup>1,2</sup> M. Bellanger,<sup>3</sup>
and L. Guillou-Frottier<sup>1,2</sup>

Received 30 March 2022; accepted 27 August 2022

Crustal fault zones provide interesting geological targets for high-temperature geothermal energy source in naturally deep-fractured basement areas. Field and laboratory studies have shown the ability of these systems to let fluid flow down to the brittle-ductile transition. However, several key questions about exploration still exist, in particular the fundamental effect of tectonic regimes on fluid flow in fractured basement domains. Based on poro-elasticity assumption, we considered an idealized 3D geometry and realistic physical properties. We examined a model with no tectonic regime (benchmark experiment) and a model with different tectonic regimes, namely a compressional, an extensional and a strike-slip tectonic regime. Compared to the benchmark experiment, the results demonstrate that different tectonic regimes cause pressure changes in the fault/basement system. The tectonic-induced pressure changes affect convective patterns, onset of convection as well as the spatial extent of thermal plumes and the intensity of temperature anomalies. Driven by poro-elastic forces, temperature anomalies around vertical faults in a strike-slip tectonic regime have a spatial extent that should be considered in preliminary exploratory phases.

**KEY WORDS:** Crustal Fault Zone, Geothermal energy, Tectonic regime, Poro-elasticity driven force, 3D dynamic THM numerical modeling.

## INTRODUCTION

Usually, the formation of a geothermal resource requires the presence of fluid, a drain and a heat source. Permeable faults are natural drains for fluids circulation and attract rising scientific and economic interest for mineral (e.g., lithium) and geothermal resources (Saevarsdottir et al., 2014; Liang et al., 2018; Guillou-Frottier et al., 2020). Crustal fault zones (CFZs) are hundreds to thousands of meter-wide geological structures that localize deformation (Ben-Zion & Rovelli, 2014) and modify the petro-physical properties of the crust from the surface down to the brittle-ductile transition (BDT). CFZs are defined by faults intersection and interconnected

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

<footnote n="2"><sup>2</sup>Univ. Orléans, CNRS, BRGM, ISTO, UMR7327, 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>Division Research/International, BASE, The Federal Office for the Safety of Nuclear Waste Management, Berlin, Germany.</footnote>

<footnote n="5"><sup>5</sup>Institute of Geological Sciences, Hydrogeology Group, Freie Universität Berlin, Berlin, Germany.</footnote>

<footnote n="6"><sup>6</sup>IFPEN, 1-4 avenue Bois-Préau, 92852 Rueil-Malmaison, France.</footnote>

<footnote n="7"><sup>7</sup>To whom correspondence should be addressed; e-mail: hugo-duwiquet@ifpen.fr</footnote>
fractures (Chester & Logan, 1986; Schulz & Evans, 1998). Field and laboratory observations show a succession of damage zones and fault cores corresponding to the multiple fault core conceptual model initially described by Faulkner et al. (2003). CFZs are present around the globe. In a few examples, they are found in Chile—Atacama Fault System (Mitchell and Faulkner, 2009), Germany–Badenweiler–Lenzkirch Suture (Brockamp et al., 2015), or Finland–Pasmajärvi Fault Zone (Ojala et al., 2019). These structures can be located in different geological settings (magmatic or amagmatic) and different tectonic stress regimes (extensional, compressional, strike-slip). This paper investigates the role of tectonic regimes on geothermal reservoirs within and around CFZs in amagmatic systems.

Although CFZs provide potential renewable and economic geothermal resources, they remain largely unexplored and therefore undeveloped. A thorough understanding of the processes that impact fluid flow and heat transfer is a prerequisite for successful geothermal exploration (Rowland & Sibson, 2004). Fluid-filled fractures have been observed down to mid-crustal depths, as shown by the deep boreholes of the Kola Peninsula (Kozlovsky, 1984), and the German KTB continental deep drill hole (Grawinkel & Stöckhert, 1997; Ito & Zoback, 2000). Additionally, Famin et al. (2004) suggested that massive infiltration of surface-derived fluids occurred in a detachment shear zone of Tinos Island (Greece), down to a depth of 10–15 km. Similarly, Siebenaller et al. (2013) demonstrated that meteoric fluid infiltration occurs down to the BDT (around a depth of 8 km) in the Naxos detachment fault (Greece). Recent work has shown that pressure and temperature conditions can induce fluid flow close to the BDT (Violay et al., 2017; Watanabe et al., 2017, 2021).

The presence of hot fluids, driven by convection around geothermal wells, is mandatory to reach production power that leverage the high drilling cost of deep borehole. When the natural permeability is insufficient, for reaching economically viable injection/production flowrates during the plant operational phase, the strategies developed by the Enhanced Geothermal System (EGS) technique consists in increasing the permeability by different injection phases, among which hydraulic stimulations aim at causing hydro-shearing and/or hydrofracturing (e.g., Gischig and Preisig, 2015; Bijay and Ghazanfari, 2021). These methods have allowed the development of many geothermal power plants,

such as Soultz-sous-Forêts, Upper Rhine Graben, France (Genter et al., 2010). However, the seismicity induced by these injection phases have jeopardized several geothermal projects (Evans et al., 2005; Deichmann & Giardini, 2009). Between 2020 and 2021 in Alsace (France) a series of induced earthquakes of magnitude *M* = 3, *M* = 3.6, *M* = 3.9 caused the definitive shutdown of the Vendenheim geothermal project and raised doubts in the population. Since 2020, the United Downs Deep Geothermal Project (UDDGP) is attempting to target a naturally fractured reservoir in the heat-producing Cornish granite (Ledingham et al., 2019; Paulillo et al., 2020). Drilling through the sub-vertical, strike-slip Porthtowan fault zone has induced an earthquake of magnitude *M* = 1.5 (www.induce.dearthquakes.org).

The complex nature of the interactions between thermal (T), hydraulic (H), mechanical (M) and chemical (C) processes affect the initial state and dynamic behavior of the geothermal reservoir under natural conditions. Theoretical and numerical modeling has been used since 1945 to understand the controlling factors of fluid circulations (Horton & Rogers, 1945; Katto & Masuoka, 1967; Horne, 1979; Forster & Smith, 1989; López & Smith, 1995; O'Sullivan et al., 2001; Magri et al., 2016; Guillou-Frottier et al., 2020). Tectonic deformation has been considered as driving forces that influence fluid flow in different geological contexts (Ord & Oliver, 1997; Cox, 1999; Rowland & Sibson, 2004). Bethke (1985) shows that compressive tectonic settings can lead to increased fluid pressure and favor upward movement. Sibson (1987) links overpressure and upward movement to fault-valve activity. In extensional tectonic settings, the generated under-pressure appears to cause downward fluid migration (McLellan et al. 2004). Cui et al. (2012) used simplified 2D models with a fault zone to show that with degrees of shortening exceeding 1% fluid flow is affected. In 3D models, Eldursi et al. (2020) suggest that during tectonically active periods, the decrease in pore pressure can reorient fluid flow in fractured zones. Nevertheless, none of them have looked at the effect of tectonic regime on fluid flow and on the temperature anomalies in CFZs.

In CFZs, fluid circulation driven by buoyancy forces occurs through upward and downward movement localizing positive temperature anomalies at shallow depths (Duwiquet et al., 2019; Guillou-Frottier et al., 2020; Duwiquet, 2022). These 2D and 3D TH (thermal-hydraulic) numerical modeling
studies have shown that vertical or subvertical deep
deformation zones could concentrate the most
important temperature anomalies at the lowest
depths. However, these results do not consider
mechanical effects on fluid flow. Among mechanical
processes, the poro-elasticity hypothesis describes
the interaction between fluids and deformation in
the porous medium. Fluids in a reservoir are af-
fected by stresses, whether on their pressure (un-
drained conditions in low-permeable media with,
e.g., increase in pressure under compressive stress
state), or on their circulation (drained conditions in
permeable media with, e.g., convection from more
to less compressed regions). Our previous study,
which included poro-elasticity, suggested that verti-
cal deformation zones oriented at 30° and 70° to a
maximum horizontal stress could correspond to
potential targets for high-temperature geothermal
energy (Duwiquet et al., 2021a; 2022). Thus, the
stress orientation has an effect on fluid flow, as al-
ready suggested by Jiang et al. (2019). In anisotropic
boundary conditions, these effects are expected to
be strongly accentuated. Indeed, they introduce
shearing conditions favorable to dilation. These as-
pects could be treated in mechanical terms, which
would allow the calculation of slip tendency, defined
by the ratio of shear stress to effective normal stress
(Morris et al., 1996).

Here, we investigate the role of tectonic re-
gimes on the formation of an amagmatic geothermal
reservoir (i.e., the only heat source is the natural
geothermal flux) within a CFZ. We propose a faul-
ted 3D THM (thermal-hydraulic-mechanical)
numerical model with a simplified geometry. The
physical properties are realistic and the permeability
is adapted to the fractured environment. In order to
understand the role of tectonic regimes, we first
considered a benchmark experiment that neglects
tectonic stress. This result is compared to several
numerical experiments where tectonic stresses are
implemented.

The effect of the poro-elasticity driven force on
fluid flow will differ from one tectonic regime to an-
other and will impact the formation of a high-temper-
ature geothermal reservoir in fractured environment.
The observed differences can be explained mainly by
different lateral fluid pressure variations. This study
highlights the importance of mechanical effects in
geothermal processes and the relevance dealing with
these aspects during the exploratory phase of such
naturally potential reservoir. The tectonic-induced
pressure changes affect fluid flow patterns, onset of

convection, spatial extent of thermal plumes and
intensity of the temperature anomalies.

METHODS

The Comsol Multiphysics™ software is based on
the finite element method (FEM) and can model,
among other various physical processes, fluid flow,
heat transfer and elastic deformation of materials in
a 3D geometry. The Comsol Multiphysics™ is a well-
known tool offering a complete access to the solu-
tion of partial differential equations. Benchmark
tests have already been performed with previous
numerical codes (OpenGeoSys, Comsol Multi-
physics™ v3.5a, and Comsol Multiphysics™ v5.4, see
Guillou-Frottier et al., (2020) and Appendix 1).

We considered an idealized and synthetic
model that represents a typical vertical fault zone of
400 m thickness, in the middle of a 5.5 km side cubic
volume of homogeneous basement rock. The fault
zone, which corresponded to a multiple fault core,
was treated like continuous porous medium. At very
high fracture density, this assumption is reasonable
(Zareidarmiyan et al., 2020). The mesh was defined
by 15,785 tetrahedra, with mesh sizes of 500 m for
low permeability zones and mesh sizes of 170 m for
high permeability zones. Preliminary convergence
tests showed that a finer mesh size gave the same
results. The transient simulations were run up to
990 kyrs. The vertical faces of the basement were
thermally insulated and fluid circulation was
blocked. On the upper horizontal face, a tempera-
ture of 20 °C and a pressure of 10⁵ Pa were imposed
(Fig. 1).

The initial thermal regime corresponded to a
geothermal gradient of 30 °C/km. At the base of the
model a heat flux of 100 Mw m⁻² was imposed. This
heat flux represents the sum of the mantle heat flow
and the heat emitted by the decay of radioactive
elements in the underlying crust. The fluid was as-
sumed to be pure water, and the fluid density de-
pended on pressure and temperature conditions (as
detailed in Duwiquet et al., 2021a). The details on
Darcy law, heat equation and Hooke law can be
found in Duwiquet et al., (2021a). All thermal, hy-
draulic and mechanical parameters are detailed in
Table 1.

The Young's moduli imposed in numerical mod-
els are suitable for fault zone and basement (Cappa &
Rutqvist, 2011). In the basement, the Poisson's ratio is
0.25, a value generally accepted in the literature.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Model setup, boundary conditions and results of the benchmark experiment (+), i.e., without stress application. Pore pressure (white contour), temperature (color patterns) and flow field (vectors) are displayed. The temperature profile is plotted against depth (right chart, crosses). When free convection is established, the 150 °C isotherm is located at 1.8 km depth. The results with stress application are shown in Figs. 4 and 5 along the ZX and YX planes. The benchmark experiment reached the steady-state regime at 65 kyrs.</figcaption></figure>

<table><caption>**Table 1.** Set of physical parameters used in numerical modeling</caption><thead><tr><th>Category</th><th>Symbols</th><th>Fault zone</th><th>Basement</th><th>Unit</th></tr></thead><tbody><tr><td>Porosity</td><td>Φ</td><td>0.1</td><td>0.05</td><td>—</td></tr><tr><td>Permeability</td><td>k</td><td>Variable</td><td>10<sup>-16</sup></td><td>m<sup>2</sup></td></tr><tr><td>Thermal conductivity</td><td>λ<sub>s</sub></td><td>3</td><td>2</td><td>W/(m.K)</td></tr><tr><td>Heat capacity</td><td>Cp<sub>s</sub></td><td>800</td><td>800</td><td>J/(kg.K)</td></tr><tr><td>Bulk density</td><td>ρ<sub>s</sub></td><td>2700</td><td>2700</td><td>kg/m<sup>3</sup></td></tr><tr><td>Young's modulus</td><td>E</td><td>5</td><td>60</td><td>GPa</td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>0.30</td><td>0.25</td><td>—</td></tr><tr><td>Biot-Willis coefficient</td><td>α<sub>B</sub></td><td>0.8</td><td>0.8</td><td>—</td></tr></tbody></table>

However, in the fault zone we imposed a value of 0.30. This value can be explained by two aspects, the first one is that the samples on which the Poisson coefficients were measured were much smaller than the block considered. This change of scale tends to increase the value of the Poisson coefficient (Heap et al., 2020). In addition, weathering and fracturing processes present in fault zones tend to increase the value of the Poisson's ratio (Heap et al., 2020). The physical properties of the fluids were identical to the study of Duwiquet et al., (2021a), as well as the incorporation of the stress evolution with depth as boundary conditions. However, as detailed in the following paragraph, here we used Andersonian assumption to consider different tectonic regimes in the 3D dynamic numerical models.

In order to investigate the influence of tectonic regimes on fluid flow, we compared a model with no tectonic stresses applied (which we called the

“benchmark experiment”) with models where different tectonic regimes were considered. For these models, we had free boundary condition at the top, and clamping at the bottom (displacement blocked in all three directions). For the four vertical sides of the model, stress boundary conditions were applied. Because we intended to study the impact of the tectonic regimes, the magnitudes of the boundary stresses were set so as to be representative of the three main tectonic regimes namely, compressional, extensional and strike-slip. We thus built up three numerical model, one per tectonic regime. For this, we considered an Andersonian assumption whereby the principal stresses are expressed with vertical (S<sub>v</sub>), maximum horizontal (S<sub>Hmax</sub>) and minimum horizontal (S<sub>hmin</sub>) components, which are used regularly in geomechanical studies of reservoirs (Anderson, 1905; Zoback et al., 2003): The relative
Tectonic Regime as a Control Factor for Crustal Fault Zone

magnitudes of these stresses determine the modeled tectonic regime, thus:

* Compressional (reverse/thrust faulting), with $S_{Hmax} \ge S_{hmin} \ge S_{V}$

* Extensional (normal faulting) with $S_{V} \ge S_{Hmax} \ge S_{hmin}$

* Strike-slip, with $S_{Hmax} \ge S_{V} \ge S_{hmin}$

Assuming that our model was aligned with the principal stresses, pure normal stresses were applied on the lateral boundaries (no shear applied on the boundaries). For the compressional regime, the fault was perpendicular to the $S_{Hmax}$ stress applied on the boundaries. For the extensional regime, the fault was perpendicular to $S_{hmin}$. Finally, for the strike-slip one, the fault was at 45° between $S_{hmin}$ and $S_{Hmax}$ orientation.

In addition, and to understand the possible effects of stress intensity, stress-depth profiles were collected on two natural systems, one in the French Massif Central and one near the San Andreas Fault. Two cases were considered, a low stress intensity zone (e.g., French Massif Central) and a high stress intensity zone (e.g., San Andreas Fault). It should be noted that we assumed that regional stresses prevail over local stress variations, which is consistent with our homogeneity assumption. The stress profiles as a function of depth corresponded to the following equations (Cornet & Burlet, 1992; Zoback, 1992):

High stress intensity:

$$ 
\begin{aligned}
\sigma_3 &= 1 \times 10^6 + (15000 \times (-z)) \\
\sigma_2 &= 4 \times 10^6 + (20000 \times (-z)) \\
\sigma_1 &= 7.5 \times 10^6 + (28667 \times (-z))
\end{aligned} 
\quad (1) $$

Low stress intensity:

$$ 
\begin{aligned}
\sigma_3 &= 0.1 \times 10^6 + (13975 \times (-z)) \\
\sigma_2 &= 0.1 \times 10^6 + (26225 \times (-z)) \\
\sigma_1 &= 0.1 \times 10^6 + (28725 \times (-z))
\end{aligned} 
\quad (2) $$

where $\sigma_1$, $\sigma_2$, $\sigma_3$, are the maximum, intermediate and minimum principal stresses (MPa), respectively; z is the vertical upwards axis, and an increase in depth (z

<0) brings a more compressive stress (positive compression convention). We assumed constant evolution of stress with depth. Notice that, in their natural state, a strike-slip stress regime holds for both regions, i.e., $\sigma_1 = S_{Hmax}$, $\sigma_2 = S_v$, $\sigma_3 = S_{Hmin}$. In the following, we kept the realistic principal stress magnitude given in Eqs. 1 and 2 but changed in our different scenarios the axes along which they operate. This way, we were able to investigate the effect of tectonic regimes while keeping realistic stress ratios. In the end, six numerical models were built (Table 2) in addition to the benchmark experiment (with no tectonic stresses applied).

Notice that in the remainder of this paper, the color coding used in Eq. 1 and is the same throughout the study for a better reading of the results. The red color code corresponds to high stress intensity and the green color code corresponds to low stress intensity. In order to allow the convergence of the 3D numerical calculations with this THM coupling, the application of the stresses was applied progressively from t = 0 yr until t = 10 yr. The results are shown in steady-state, at 65 kyr for the benchmark experiment and at 990 kyr for the model where stresses were applied.

Fluid flow velocities, and thus efficiency of heat transport by convection through the crust, were controlled directly by the permeability of the rocks. Permeability is a dynamic and variable parameter that can change during different geological processes (Gleeson & Ingebritsen, 2016). The relationship between permeability and porosity is widely debated. However, although adapted to porous media, the Kozeny-Carman relationship does not seem to be appropriate for fractured media (Lamur et al., 2017; Parisio et al., 2019). Instead, the permeability k we considered has been empirically demonstrated to be appropriate for fractured media (Lamur et al., 2017), thus:

$$ 
\begin{aligned}
\log k &= (1 - \omega)\log k_i + \omega\log k_f \\
k_i &= 4.979 \times 10^{-11} n^{3.11} \\
k_f &= 1.143 \times 10^{-11} n^{0.64}
\end{aligned} 
\quad (3) $$

where $\omega$ represents the degree of rock fracturing; $\omega = 1$ for fully fractured rock and $\omega = 0$ for fully intact rock. Here, we took $\omega = 0.8$ corresponding to the significant fracturing degree of the CFZ, interspersed with intact zones. The n represents initial porosity. Here, n = 0.05 for intact rock and n = 0.1
<table><caption>**Table 2.** Synthesis of the different tectonic regimes and stress intensities considered</caption><thead><tr><th>Stress intensity</th><th>Tectonic regime</th><th>Normal stress applied on boundaries</th></tr></thead><tbody><tr><td rowspan="3">High (Eq 1)</td><td>Extensional</td><td>Boundaries parallel to fault orientation : S<sub>hmin</sub><br>Boundaries perpendicular to fault orientation : S<sub>Hmax</sub></td></tr><tr><td>Compressional</td><td>Boundaries parallel to fault orientation : S<sub>Hmax</sub><br>Boundaries perpendicular to fault orientation : S<sub>hmin</sub></td></tr><tr><td>Strike-slip</td><td>Fault orientation oblique towards lateral boundaries</td></tr><tr><td rowspan="3">Low (Eq 2)</td><td>Extensional</td><td>Boundaries parallel to fault orientation : S<sub>hmin</sub><br>Boundaries perpendicular to fault orientation : S<sub>Hmax</sub></td></tr><tr><td>Compressional</td><td>Boundaries parallel to fault orientation : S<sub>Hmax</sub><br>Boundaries perpendicular to fault orientation : S<sub>hmin</sub></td></tr><tr><td>Strike-slip</td><td>Fault orientation oblique towards lateral boundaries</td></tr></tbody></table>

for fractured rock. The $k_i$ is permeability (m²) of intact rock, and $k_f$ is permeability (m²) of fractured rock. The stress state of a fault can be defined qualitatively in terms of slip tendency, defined by the ratio of shear stress to effective normal stress (Morris et al., 1996). In purely elastic domain and without dissipative phenomena, the slip tendency is an indicator (high/low general trend) to assess the stress variation on a fault zone. More specifically, slip tendency indicates how critically stressed fault zones are, and how easily they can be reactivated by, e.g., stimulation methods. Consequently, some preferential flow paths may appear under the action of a given applied tectonic stress (Siler et al., 2019). The slip tendency is defined as:

$$T_s = \frac{\tau}{\sigma'_n} \qquad (4)$$

where $\tau$ is shear stress and $\sigma'_n$ is effective normal stress defined as:

$$\sigma'_n = \sigma_n - (\alpha_B \times pf) \qquad (5)$$

where $\sigma_n$ is total normal stress; $\alpha_B$ is the Biot-Willis coefficient and $pf$ is fluid pressure.

# RESULTS

## Formation of a High-Temperature Geothermal Reservoir Within a Vertical Crustal Fault Zone with No Tectonic Stresses: The Benchmark Experiment

The numerical simulation without stress application showed an upward circulation at the fault

center and two downward circulations at the fault ends (Fig. 1). As illustrated by colored streamlines in Figure 1, these fluid circulations transferred heat by convection. Thermally induced buoyancy forces drove the flow. The hot and less dense fluid rose to the surface. As its density increases due to the cooling effects of the top boundary condition, the fluid sank within the fault plane leading to a classical natural convection pattern. Under these conditions, the free convection generated a thermal disturbance of the environment. In a purely diffusive setting, the 150 °C isotherm was observed at 5 km depth, whereas here the same isotherm was localized at 1.8 km depth, at the center of the fault. Fluid velocities were of the order of $1 \times 10^{-9}$ m.s$^{-1}$ at the center of the fault and $20 \times 10^{-9}$ m.s$^{-1}$ at the ends of the fault. The higher fluid velocities at the bottom center of the fault were responsible for slight increase in pressure (white lines). This trend was already observed in other numerical TH modeling (Scott et al., 2017). The temperature profile at the center of the model was similar to typical temperature profiles measured in geothermal reservoirs such as Soultz-sous-Forêts (Genter et al., 2010).

## Application of Tectonic Stresses

The effects of high and low stress intensity application (Eqs. 1 and 2 “method section”) on the convective regimes are shown in Figures 2 and 3, respectively. The results are shown on vertical cross sections along the fault plane and on horizontal cross sections, at 2 km depth. For each tectonic regime that was considered, the results showed the temperature anomalies, the fluid flow pattern and
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Results of numerical modeling after high stress intensity stresses application (in red, see Eq. 2). The results are shown in vertical section (side view, located in the middle of the fault) and in horizontal section (top view, located at 2 km depth). The scale of temperature anomalies and fluid flow velocities is different according to the tectonic regimes. For each regime, the maximum and minimum values of temperatures and fluid flow velocities are indicated. Positive temperature anomalies are colored red, negative temperature anomalies are colored blue. Fluid circulation is marked by the lines, the direction by the arrows. The color of the lines corresponds to fluid velocity. In red, the velocity is the fastest; in blue, the velocity is the slowest.</figcaption></figure>

the fluid velocities. The results can be compared directly to the benchmark experiment (Fig. 1). Notable differences existed and are detailed below.

*Temperature Anomalies*

Regardless of the considered tectonic regime, different positive and negative temperature anomalies were observed. In the following, temperature anomalies correspond to the difference of temperature resulting from a conductive thermal regime.

They differed from the benchmark case by their number, intensity and lateral extension.

In an extensional tectonic regime, two positive temperature anomalies formed (Figs. 2 and 3). They were + 55 °C in high stress intensity and + 47 °C in low stress intensity. For both stress intensities, the second temperature anomalies were less intense, + 10 °C in high stress intensity and + 7 °C in low stress intensity. For both stress intensities, a negative temperature anomaly developed at the center of the fault. This anomaly was - 60 °C for the high stress
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Results of numerical modeling after low stress intensity application (in green, see Eq. 1). The results are shown in vertical section (side view, located in the middle of the fault) and in horizontal section (top view, located at 2 km depth). The scale of temperature anomalies and fluid flow velocities is different according to the tectonic regimes. For each regime, the maximum and minimum values of temperatures and fluid flow velocities are indicated. Positive temperature anomalies are colored red, negative temperature anomalies are colored blue. Fluid circulation is marked by the lines, the direction by the arrows. The color of the lines corresponds to fluid velocity. In red, the velocity is the fastest; in blue, the velocity is the slowest.</figcaption></figure>

intensity and – 57 °C for the low stress intensity. For a high stress intensity, at 2 km depth, the horizontal cross section showed a negative temperature anomaly that reached a maximum of – 30 °C. This anomaly covered a large surface area of the fault. The positive temperature anomaly of + 45 °C occupied the remaining space, but extended further into the basement. In the basement and up to the edge of the model, we found a positive temperature anomaly of + 20 °C. For a low stress intensity, at 2 km depth, two positive temperature anomalies and one negative anomaly were observed. The first positive anomaly was + 41 °C, the second had a

value of + 20 °C. This + 20 °C anomaly spread more widely in the system while the + 41 °C anomaly remained localized on the fault. The negative anomaly located at the center of the fault was – 37 °C and extended into the basement in the same way as the positive + 41 °C anomaly.

In a compressional tectonic regime, independently from the stress intensity, two positive temperature anomalies were observed. The value of the maximum temperature anomaly in high stress intensity was + 60 °C and in low intensity it was + 51 °C. Depending on the stress intensity, these two maximum values were spatially slightly shifted
Tectonic Regime as a Control Factor for Crustal Fault Zone

(Figs. 2 and 3). The amplitudes of the second tem-
perature anomaly were + 35 °C and + 40 °C for low
and high stress intensity, respectively. In the hori-
zontal cross section, the two temperature anomalies
were found for each intensity. In high stress inten-
sity, the values of these two temperature anomalies
were + 56 °C and + 47 °C. In low stress intensity,
the values were + 47 °C and + 35 °C. The lateral
extension of these temperature anomalies was lim-
ited. They were surrounded by negative temperature
anomalies that locally reached – 47 °C in high
intensity and – 42 °C in low intensity. At a depth of
2 km, the positive temperature anomalies were
much less extended than in extensional tectonic or
strike-slip regimes (see below).

In strike-slip regimes and regardless of stress
intensity applied (Figs. 2 and 3), positive tempera-
ture anomaly extended widely along the length of
the fault, from the surface to 4.5 km deep. The
maximum positive temperature anomaly value
was + 70 °C for the high stress intensity (Fig. 2),
and + 62 °C for the low stress intensity (Fig. 3). In
horizontal cross section, these temperature anoma-
lies spread largely beyond the fault.

Temperature anomalies of + 25 °C were found
in the basement, suggesting that in a strike-slip sys-
tem, the positive temperature anomaly, and inde-
pendently from the stress intensity, represents an
important volume. This heat propagation was done
by thermal diffusion from the fault center, where the
temperature anomaly was the most intense. Indeed,
the larger the convection cell inside the fault, the
wider the extent of the diffusive perturbation. At
2 km depth, the maximum value of the temperature
anomaly was + 65 °C for the high stress intensity,
and + 61 °C for the low stress intensity. Negative
temperature anomalies were present and localized at
the extremities of the fault. They were – 45 °C for
the high stress intensity and – 40 °C for the low
stress intensity.

To summarize, tectonic regimes influence the
distribution and the amplitude of temperature
anomalies. Positive temperature anomalies were
most intense in strike-slip, then in compression and
extension. The spatial extent of positive temperature
anomalies was not identical for each tectonic regime.
In strike-slip, these anomalies were largely extended
through the basement. This lateral extension was
less important in the extensive tectonic regime. Fi-
nally, in compressional regime these anomalies were
localized in the near vicinity of the fault. The tec-
tonic regimes have shown to play a key role in

temperature distribution, and this was clearly re-
lated to the different convective patterns and fluid
flow velocity, as described below.

Fluid Flow Pattern

In extensional tectonic regimes, the fluid flow
pattern was characterized by a downward movement
at the center of the fault and two upward move-
ments at the ends of the fault (Figs. 2 and 3). In high
stress intensity the maximum fluid velocity was
30 × 10⁻⁹ m.s⁻¹, against 21 × 10⁻⁹ m.s⁻¹ in low
stress intensity. The minimum fluid velocities were,
for each intensity, 1 × 10⁻⁹ m.s⁻¹. In horizontal
cross section, at 2 km depth, upward movements
were generated at positive temperature anomalies
while downward movements were generated at
negative temperature anomalies. In high stress
intensity, the maximum fluid velocity was 25 × 10⁻
⁹ m.s⁻¹, compared to 17 × 10⁻⁹ m.s⁻¹ in low stress
intensity. The minimum fluid velocity was 5 × 10⁻
⁹ m.s⁻¹ in high stress intensity, compared to 1 × 10⁻
⁹ m.s⁻¹ in low stress intensity.

In compressional tectonic regimes, for each
stress intensity, there was a slightly different con-
nective pattern. In high stress intensity, there were
two upward and two downward movements,
whereas in low stress intensity there were two up-
ward and three downward movements. The maxi-
mum fluid velocity for high stress intensity was
5 × 10⁻⁹ m.s⁻¹ and for low stress intensity was
10 × 10⁻⁹ m.s⁻¹. For this tectonic regime, the max-
imum fluid velocity was present when the lowest
stress intensity was applied. In horizontal cross sec-
tion, at 2 km depth and for high stress intensity,
upward movements were localized where tempera-
ture anomalies were positive, whereas downward
movements were localized where temperature
anomalies were negative. The fluid velocity was
between 1 and 2 × 10⁻⁹ m.s⁻¹. In low stress inten-
sity, the downward movements were localized at the
extremities and at the center of the fault. The up-
ward movements were localized between each
downward movement. The maximum fluid velocity
was 10 × 10⁻⁹ m.s⁻¹ and the minimum was 1 × 10⁻
⁹ m.s⁻¹. At a depth of 2 km, fluid velocities varied
from 1 to 2 × 10⁻⁹ m.s⁻¹. Thus, fluid velocity values
in a compressional regime were the lowest recorded.

In the strike-slip tectonic regime, as for the
benchmark experiment, the fluid flow pattern was
characterized by an upward movement at the center
of the fault and two downward movements at the ends of the fault (Figs. 2 and 3). The fluid velocity varied with stress intensity. In low stress intensity, the minimum fluid velocity was $1 \times 10^{-9}$ m.s$^{-1}$ and the maximum was $12 \times 10^{-9}$ m.s$^{-1}$ at the downward movements. In high stress intensity, the minimum fluid velocity was $1 \times 10^{-9}$ m.s$^{-1}$ and the maximum was $17 \times 10^{-9}$ m.s$^{-1}$.

With application of tectonic stresses, the fluid flow pattern was different from the benchmark experiment, with consequences on temperature anomalies and fluid velocities. In the benchmark experiment, buoyancy was the only driving force for fluid convection. Here, in the presence of tectonic stresses, stress induced forces also influenced fluid flow.

## Lateral Fluid Pressure Variation

After stress application, the fluid pressure was different from the benchmark experiment where pressure was hydrostatic (Fig. 4). Whatever the stress intensity applied, the experiment with compressional and strike-slip tectonic regimes showed fluid pressures higher than those in the benchmark experiment, whereas the extensional tectonic regime displayed fluid pressures lower than those in the

benchmark experiment (Fig. 4). The fluid pressure varied between the basement and the fault. This lateral variation differed according to the stress intensity applied.

In each of the cases, the tectonic regimes and the stress intensity generated lateral fluid pressure variation between the fault and the basement. In high stress intensity, and for the same horizontal distance, these lateral pressure differences were 1.5, 0.45, and 0.97 MPa for compressional, strike-slip, and extensional tectonic regimes, respectively. In low stress intensity, the lateral pressure differences were 1.35, 0.3 and 0.91 MPa for compressional, strike-slip, and extensional tectonic regimes, respectively. These lateral fluid pressure differences drove the fluids from the high-pressure zones (i.e., basement) to the low-pressure zone (i.e., fault).

## What effect(s) on the Onset of Positive Temperature?

The regional mechanical stresses imposed a pressure distribution that interacted with buoyancy forces, which caused free convection patterns. The lateral fluid pressure differences brought the fluid from the zones of high pressure to those of low pressure. However, as shown in Figure 5, the time needed to set up this temperature anomaly varied.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Lateral fluid pressure in high stress intensity application (in red) and low stress intensity application (in green), at 2 km depth. For comparison, the fluid pressure of the benchmark experiment is indicated (+). With stresses application, the fluid pressure is different from the benchmark experiment. The lateral fluid pressure variation depends on the tectonic regime and stress intensity.</figcaption></figure>
Tectonic Regime as a Control Factor for Crustal Fault Zone

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Time required to reach the steady-state temperature anomalies, depending on the tectonic regimes and the applied stress intensities. The color code and patterns are the same as those used in the previous figures.</figcaption></figure>

Considering the compressional tectonic regime (solid line), the time needed to reach a steady-state regime in high stress intensity was 85 kyrs, while it was 95 kyrs in low stress intensity. This difference was observed for each regime so that the stress intensity had a role in the onset of convection. As can be seen in Figure 5, the compressional regime, which had the highest lateral fluid pressure difference, reached the steady-state more quickly than the other tectonic regimes. The lowest lateral fluid pressure difference was recorded for the strike-slip regime. The strike-slip regime reached the steady-state regime at 203 kyrs for high stress intensity and 223 kyrs for low stress intensity.

The benchmark experiment reached the steady-state regime as early as 65 kyrs (Fig. 5), and finally the additional stress input may increase the time to reach a steady-state in situ. The differences observed on the different tectonic regimes were consistent with the differences observed on the lateral fluid pressure between the fault and the basement. However, pressure differences alone could not explain the different fluid velocities observed (Figs. 2 and 3), and the mechanical behavior of the basement-fault system after stress application. In a geothermal exploration context, these numerical calculations showed that, in a strike-slip tectonic

regime, the temperature anomalies should be more intense than in a compressive or extensive regime, but that the onset was the most delayed with respect to the benchmark case.

*Permeability and Slip Tendency Variation*

In the fault, very small permeability variations were observed as a function of stress intensity and tectonic regime. For the benchmark experiment, the fault's permeability was $1 \times 10^{-12.1} \text{ m}^2$. In high stress intensity, the variation of the permeability values in the center of the fault was as follows: 0.08%, 0.42% and 0.66% for extensional, strike-slip and compressional tectonic regimes, respectively. In low stress intensity, the variation of the permeability values in the center of the fault was as follows: 0.08%, 0.33%, 0.74% for extensional, strike-slip and compressional tectonic regimes, respectively. The effect of stress intensity on permeability remained negligible. For example, considering the compressional tectonic regime, the difference between high and low stress intensity was 2.3%. We noted a consistency between the values of fluid velocities and permeability values. Indeed, the fluid velocities were the most important for the highest permeability, and
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Slip tendency in high stress intensity application (in red), and low stress intensity application (in green).</figcaption></figure>

the least important for the lowest permeability. Once the stresses were applied, the permeability was slightly different from the benchmark experiment, and controlled the fluid velocity.

The general trend of slip tendency is shown for each tectonic regime and stress intensity applied (Fig. 6). Overall, the general trend was higher in high stress intensity than in low stress intensity applied. The results showed a distinctive variation in the fault and in the basement for each tectonic regime. Regardless of the stress intensity applied, the slip tendency in a strike-slip tectonic regime was the highest, because it was more critically stressed. In the strike-slip regime, the slip tendency increased slightly from the edge of the basement to the level of the fault zone. At the edge of the fault, the slip tendency increased significantly. In the middle of the fault, the slip tendency decreased slightly. This evolution was the same regardless of the stress intensity applied.

In compressional and extensional regimes, the trend from the edges of the basement to the fault increased similarly. This increase was slightly greater in compression than in extension. At the edge of the fault, the general trend was the same. In compression, the slip tendency decreased, and in extension it increased. In extension, the middle of the fault was less critically stressed than in compression. Regardless of the intensity of the stress, this evolution was the same. Of the three tectonic regimes tested, the strike-slip regime was the most critically stressed. The general trend showed that it is at the interface

between the basement and the fault that the trend was most important. This result illustrates the heterogeneity of the mechanical parameters incorporated in the numerical model, and it shows further that the role of the stress intensity on the range of values tested will not drastically change the mechanical response of the basement-fault system.

## DISCUSSION

The presence of fluids and a sufficiently high permeability are two factors that allow the development of a geothermal resource, from a natural geothermal flux. Considering a simplified geometry, but realistic physical properties (pressure and temperature-dependent fluid density, temperature-dependent fluid viscosity), our 3D TH and 3D THM numerical modeling showed that, in a basement domain, CFZs allow hot fluids to ascend at economically viable depth. The comparison of numerical models with and without tectonic stresses highlights the non-negligible role of tectonic regime on the spatial distribution of positive temperature anomalies.

### Different Tectonic Regime, Different Fluid Flow

Without stress application, the benchmark experiment showed an upward movement at the center of the fault and brought about a temperature
# Tectonic Regime as a Control Factor for Crustal Fault Zone

of 150 °C at 1.8 km depth (Fig. 1). Two downward movements were located at the ends of the fault. This convective pattern was already observed in 2D and 3D TH numerical modeling of fault zones in the crustal domain (Wanner et al., 2019; Guillou-Frottier et al., 2020) and it was called “bulb-like” convective pattern. This particular shape is favored because hot fluids have less resistance when flowing upwards (viscosity is smaller at high temperature) and it can also localize temperature anomalies in the basement. With a flat topography and without metamorphic devolatilization and/or magmatic fluid production that might change fluid pressure (Nur & Walder, 1990), fluid circulation is driven only by buoyancy force. With tectonic regimes application, we saw that other forces are added and modify the fluid circulation.

Considering the experiments where tectonic regimes are accounted for, the obtained results were different from the benchmark experiment (Figs. 2 and 3). This highlights the key role of the relationship between tectonic settings and fluid flow. In our results, the stress intensity parameter did not change the general dynamics of convective patterns between two identical regimes. Considering a sensitivity study on a large Pontgibaud Crustal Fault Zone (French Massif Central), where the stresses boundary magnitudes were tested while holding their directions, a shift from one convective pattern to another would occur for a variation of approximately 40 MPa in the maximum horizontal stress magnitude (Duwiquet et al., 2021a; 2021b). For this last example and in the light of the observations made in this study, it could be possible that another force besides buoyancy influences fluid circulation. In this study, the stress intensity applied was not a factor influencing the convection pattern. However, the stress intensity has a role in the time to reach the steady state. Indeed, the steady-state was reached faster when the most intense stresses were applied (Fig. 5), and this was not explored in the Duwiquet et al. (2021a; 2021b) study. Moreover, the role of tectonics on the direction of fluid flow has been highlighted. Considering a numerical approach, studies have shown that, in compression, upward fluid movements are favored (Upton, 1998), whereas, in extension, downward fluid movements are favored (Cui et al., 2012). Here, we saw that different tectonic regimes have a role in the onset of convection.

# Nature and Impact of Poro-elasticity Driven Force

These numerical models showed that the tested tectonic regimes had an influence on fluid pressure variation. The poro-elastic assumption describes the interaction of fluids and deformation in porous media. Compressional, extensional, and strike-slip tectonic regimes induce variable fluid pressures for each case. The incorporation of heterogeneous mechanical parameters between the fault and the basement (see Table 1) led to a different mechanical response and thus to a heterogeneous variation of the fluid pressure between the fault and the basement (see Fig. 4). This lateral fluid pressure difference drove the fluids from the high-pressure zones to the low-pressure zones. This effect was facilitated by higher permeability values in the fault than in the basement.

In compressional tectonic regime, the positive temperature anomaly was concentrated in the fault, whereas in strike-slip tectonic regime, this anomaly extended more widely into the basement (Fig. 4). This is related to the different fluid pressures. The difference in fluid pressure between the fault and the basement depended on the tectonic regime. In compressional tectonic regime, the difference in fluid pressure between the basement and the fault was larger than in strike-slip tectonic regime. This concentrated the fluid flow in a more restricted volume in compression than in strike-slip. This impacts the spatial extent of the positive temperature anomaly, and therefore the geothermal potential of these naturally fractured zones.

The force that drives fluids from high-pressure zones to low-pressure zones is similar to the effects of topography on fluid flow (Forster & Smith, 1989; López & Smith, 1995). This poro-elasticity driven force will therefore influence the initial fluid motions (downward and upward movement) and thus the time required for the steady-state temperature anomaly to develop. Considering a TH coupling, the fluid circulation in the benchmark experiment was driven by free convection. Considering three different tectonic regimes with the same TH coupling, the fluid circulation was the result of an interplay between forced and free convection. These changes were not related to boundary conditions, but to the tectonic regimes themselves. Forced convection is here referred to as stress induced convection.

After stress applications, the permeability values obtained were consistent with fractured and altered granitic environment (Sardini et al., 1997;
Duwiquet et al., 2019; 2021a; Gomila et al., 2021). This should be sufficient to allow fluids to flow by buoyancy forces and transfer heat by free convection. However, between tectonic regimes, the difference in permeability cannot explain convective pattern variability results. Nevertheless, these small permeability variations influence the velocities of fluid flow. The fluid velocity decreased with permeability decrease (see Figs. 2, 3, and sub-section *Permeability and Slip Tendency Variation*). Actually, the positive temperature anomalies in the compression and extensional regime model remained centered on the fault compared to the strike-slip tectonic regimes. Therefore, other processes must limit the development of temperature anomalies in compressional regime. By concentrating the flow of fluids on limited spaces, the poro-elasticity driven force could have the effect of spatially concentrating temperature anomalies when the difference between the fluid pressure of the basement and the fault is large.

## Consequences for High-Temperature Geothermal Potential of Crustal Fault Zones

For geothermal exploration, slip tendency analyses can be used to target favorable zones for natural fluid flow and future enhancement by fault reactivation (Barton et al., 1995; Morris et al., 1996; Ito & Zoback, 2000). Our results showed that the strike-slip regime would be the most favorable to allow fluid flow. The occurrence of geothermal reservoirs in such contexts was already known, as for the Alpine Fault, New Zealand (Boulton et al., 2012), or in the Geysers geothermal area, California (Altmann et al., 2013). In the light of our results, an exploratory phase on the geothermal potential of fault zones could further consider strike-slip faulted regime as interesting targets for geothermal power system.

If slip tendency can be used as a potential qualitative indicator in purely elastic models, a more accurate interpretation would require incorporating dissipative mechanical behaviors, such as e.g., Mohr-Coulomb elasto-plastic law, which would allow to quantify the variation in slip tendency. By considering irreversible mechanical processes, such as including dilation (opening under shearing) and fracturing (increase in pore space/fractures), the highest favorability of strike-slip regime toward convection is expected to be emphasized.

We used an idealized geometry with single, vertical fault zone. Hence, the linear stress boundary conditions give rather homogeneous stress ratios and states along the fault. In natural systems, the network of variously oriented and dipping fault zones brings heterogeneous stress states, even along each fault zone taken separately. Still, and recalling that the aim of our study is to better understand the setup conditions of convection cells, the results already allow us to highlight the complex impact of mechanics over convection patterns. Among others, our conclusions tend to show that within a complex fault zone network, fault zones undergoing strike-slip conditions might be the ones to be preferably explored. Future work, including sensitivity analysis on parameters, law of behavior, boundary condition effects, should be conducted to confirm this trend.

These fundamental results are generally applicable in nature to any fracture rock that may contain fluids, gas or oil in its pores, in basement geological context. However, our results do not take into account fundamental aspects such as the consideration of more complex rheology, to start with plastic phenomena. The effects of fault intersections, precipitation and dissolution of mineral phases, and fluid composition were also not taken into account. Moreover, temperature impact on stress distribution around a fault should be investigated. These interrelations induce shear stresses that can reactivate a fault (Karaoğlu et al., 2020), and thereafter modify its permeability. The dependence of permeability on all of these phenomena should probably generate permeability anisotropy. The development of this anisotropy can cause a change in the intrinsic properties of a geothermal reservoir and induce a change in the heat transfer mode (Sun et al., 2017). In cases where fluid flow is important (i.e., areas of high permeability), the permeability can take a particular geometry that minimizes the resistance to flow and then optimizes fluid flow (Bejan & Lorente, 2011). Fluid salinity could also play an important role, but effects would be marked above 400 °C (Driesner, 2007), which was beyond our modeled temperatures.

Geothermal energy can become a major asset for the transition to low-carbon energy sources. Its development requires, prior to exploration, comprehensive understanding of limiting and enabling factors controlling fluid flow and the location of temperature anomalies (Jolie et al., 2021). Fault intersections, triple junctions systems, are, among others, examples of anomalously permeable zones,
Tectonic Regime as a Control Factor for Crustal Fault Zone

that localize fluids flow (Person et al. 2012; Karaoğlu et al., 2016, 2019). Our numerical results showed that, within a single permeable fault, and without other heat sources, the poro-elasticity driven force provoked by the tectonic regimes causes lateral fluid pressure variation, allowing for the more or less rapid development of the temperature anomaly, by mixed, free and forced convection. For example, strike-slip tectonic regimes would have the largest temperature anomalies, but would take the longest time to set up.

Overall, this work suggests that anomalous permeable zones, like CFZs, with no external heat source, have significant energy potential. The exploration of vertical CFZs in strike-slip stress regime could accelerate the transition to low-carbon, renewable and climate-neutral energy. Studies have shown that major strike-slip fault zones localize porosity and permeability even beyond the BDT (Faulkner et al., 2010; Cao & Neubauer, 2016). Our study confirms that this tectonic regime seems to favor higher thermal anomalies than compressional tectonics.

The complex nature of the processes that occur during the development of a geothermal resource within a CFZ makes it an environment that requires a state-of-the-art numerical analysis of the processes that can control fluid circulation. Such fully coupled analysis will help exploration phases and ultimately promote the development of this renewable energy, not only in anomalously hot areas, but in anomalous permeable areas, and thus promote the development of this climate-neutral energy.

FUNDINGS

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), TLS-Geothermics and by the "GERESFAULT" project (ANR-19-CE05-0043).

ACKNOWLEDGMENTS

HG warmly thanks the DGR/GEM team of BRGM for their hospitality during the preparation of this manuscript as well as Institut Français du Pétrole et des Energies Nouvelles (IFPEN) for permission to submit this manuscript for publication.

HD would also like to sincerely thank Patrick Ledru, Frédéric Donzé, Roland Horne, Christine Souque, Albert Genter and Fabrice Gaillard who, through their relevant advice and suggestions within the framework of the thesis defense, have contributed to improve the final version of this manuscript. All authors would like to thank the anonymous reviewer (R1) as well as Özgür Karaoğlu (R2) for their insightful suggestions as well as their constructive remarks

DECLARATIONS

**Conflict of Interest** The authors declare that they have no conflict of interest.

OPEN ACCESS

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

APPENDIX 1

The calibration of our numerical experiments was performed based on results from Comsol Multiphysics™ version 3.5 (Guillou-Frottier et al., 2020) as well as the OpenGeoSys numerical code (Magri et al., 2017). These models consider a 40 m wide vertical fault in an impermeable box (a 5.5 × 5.5 × 5.5 km side cube). The fluid properties were identical for all three results with a linear dependence of temperature with water density, and an exponential decrease in viscosity with temperature. This result is shown for a time $t_0 + 10^{13}$ s. The imposed permeability value was $5 \times 10^{-15} m^2$. The
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Benchmarking of our numerical experiment with the OpenGeoSys Code (Magri et al., 2017) and Comsol Multiphysics.™ V3.5a (Guillou-Frottier et al., 2020).</figcaption></figure>

fluid flow velocity was $2.03 \times 10^{-9}$ m.s$^{-1}$. The fluid flow velocity was slightly higher than that described by Guillou-Frottier et al., (2020), who recorded a velocity of $1.4 \times 10^{-9}$ m.s$^{-1}$.

The convective patterns were similar. Fluid velocity accelerated along the permeable fault and exhibited upward movement due to the thermal gradient and buoyancy forces related to lower water density at depth. The thermal perturbations were also within the orders of magnitude of previous studies with - 23.04 °C and + 31.14 °C.

The numerical experiments of Magri et al. (2017) and Guillou-Frottier et al. (2020) used Darcy's law in conjunction with the equations of heat. The numerical experiments in this study coupled the equations of heat and mass transfer with the modulus of poro-elasticity. The poro-elasticity interface of the Comsol Multiphysics™ combines Darcy's law with the linear elastic behavior of porous media (details can be found in Duwiquet et al. 2021a). Poro-elastic coupling allows boundary stresses to be imposed, which can be recorded as fluid pressure. In Figure 7 (right panel), no stress or other mechanical conditions were applied. Consequently, this experiment corresponded to an identical coupling as in the previous studies (Fig. 7, left and middle panels).

## REFERENCES

Altmann, J. B., Heidbach, O., & Gritto, R. (2013). Relative importance of processes leading to stress changes in the

Geysers geothermal area. In *Proceedings of Thirty-Eighth Workshop on Geothermal Reservoir Engineering*.

Anderson, E. M. (1905). The dynamics of faulting. *Transactions of the Edinburgh Geological Society*, 8(3), 387–402.

Barton, C. A., Zoback, M. D., & Moos, D. (1995). Fluid flow along potentially active faults in crystalline rock. *Geology*, 23(8), 683–686.

Bejan, A., & Lorente, S. (2011). The constructal law and the evolution of design in nature. *Physics of Life Reviews*, 8(3), 209–240.

Ben-Zion, Y., & Rovelli, A. (2014). Properties and processes of crustal fault zones: volume I. *Pure and Applied Geophysics*, 171(11), 2863–2865.

Bethke, C. M. (1985). A numerical model of compaction-driven groundwater flow and heat transfer and its application to the paleohydrology of intracratonic sedimentary basins. *Journal of Geophysical Research: Solid Earth*, 90(B8), 6817–6828.

Bijay, K. C., & Ghazanfari, E. (2021). Geothermal reservoir stimulation through hydro-shearing: An experimental study under conditions close to enhanced geothermal systems. *Geothermics*, 96, 102200.

Boulton, C., Carpenter, B. M., Toy, V., & Marone, C. (2012). Physical properties of surface outcrop cataclastic fault rocks, Alpine Fault, New Zealand. *Geochemistry, Geophysics, Geosystems*. https://doi.org/10.1029/2011GC003872.

Brockamp, O., Schlegel, A., & Wemmer, K. (2015). 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*, 104(3), 683–702.

Cao, S., & Neubauer, F. (2016). Deep crustal expressions of exhumed strike-slip fault systems: Shear zone initiation on rheological boundaries. *Earth-Science Reviews*, 162, 155–176.

Cappa, F., & Rutqvist, J. (2011). Impact of CO2 geological sequestration on the nucleation of earthquakes. *Geophysical Research Letters*. https://doi.org/10.1029/2011GL048487.

Chester, F. M., & Logan, J. M. (1986). Implications for mechanical properties of brittle faults from observations of the Punchbowl fault zone, California. *Pure and applied geophysics*, 124(1), 79–106.

Cornet, F. H., & Burlet, D. (1992). Stress field determinations in France by hydraulic tests in boreholes. *Journal of Geophysical Research: Solid Earth*, 97(B8), 11829–11849.
Tectonic Regime as a Control Factor for Crustal Fault Zone

Cox, S. F. (1999). Deformational controls on the dynamics of fluid flow in mesothermal gold systems. *Geological Society, London, Special Publications, 155*(1), 123-140.

Cui, T., Yang, J., & Samson, I. M. (2012). Tectonic deformation and fluid flow: Implications for the formation of unconformity-related uranium deposits. *Economic Geology, 107*(1), 147-163.

Deichmann, N., & Giardini, D. (2009). Earthquakes induced by the stimulation of an enhanced geothermal system below Basel (Switzerland). *Seismological Research Letters, 80*(5), 784-798.

Driesner, T. (2007). The system H2O-NaCl. Part II: Correlations for molar volume, enthalpy, and isobaric heat capacity from 0 to 1000 C 1 to 5000 bar, and 0 to 1 XNaCl. *Geochimica et Cosmochimica Acta, 71*(20), 4902-4919.

Duwiquet, H. et al. (2021b). Crustal fault zones (CFZ) as geothermal power systems: 3D variation of permeability and related processes. In *Stanford Geothermal Workshop*, (2021b, February).

Duwiquet, H. (2022). Crustal fault zones as geothermal power systems. *Contribution of Numerical Modelling and Comparison with Natural Systems* (Doctoral dissertation, Université d'Orléans).

Duwiquet, H., et al. (2022). Crustal Fault Zones as underexploited geothermal resources: Contribution of numerical modelling and comparison with natural systems. In *Stanford Geothermal Workshop*, (2022, February).

Duwiquet, H., et al. (2019). On the geothermal potential of crustal fault zones: A case study from the Pontgibaud area (French Massif Central, France). *Geothermal Energy, 7*(1), 1-29.

Duwiquet, H., et al. (2021a). Crustal fault zones (CFZ) as geothermal power systems: A preliminary 3D THM model constrained by a multidisciplinary approach. *Geofluids*. https://doi.org/10.1155/2021/8855632.

Eldursi, K., et al. (2020). New insights from 2-and 3-D numerical modelling on fluid flow mechanisms and geological factors responsible for the formation of the world-class Cigar Lake uranium deposit, eastern Athabasca Basin, Canada. *Mineralium Deposita, 56*(7), 1365-1388.

Evans, K. F., et al. (2005). Microseismicity and permeability enhancement of hydrogeologic structures during massive fluid injections into granite at 3 km depth at the Soultz HDR site. *Geophysical Journal International, 160*(1), 388-412.

Famin, V., Philippot, P., Jolivet, L., & Agard, P. (2004). Evolution of hydrothermal regime along a crustal shear zone, Tinos Island, Greece. *Tectonics*. https://doi.org/10.1029/2003TC001509.

Faulkner, D. R., et al. (2010). A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones. *Journal of Structural Geology, 32*(11), 1557-1575.

Faulkner, D. R., Lewis, A. C., & Rutter, E. H. (2003). On the internal structure and mechanics of large strike-slip fault zones: Field observations of the Carboneras fault in south-eastern Spain. *Tectonophysics, 367*(3-4), 235-251.

Forster, C., & Smith, L. (1989). The influence of groundwater flow on thermal regimes in mountainous terrain: A model study. *Journal of Geophysical Research: Solid Earth, 94*(B7), 9439-9451.

Genter, A., Evans, K., Cuenot, N., Fritsch, D., & Sanjuan, B. (2010). Contribution of the exploration of deep crystalline fractured reservoir of Soultz to the knowledge of enhanced geothermal systems (EGS). *Comptes Rendus Geoscience, 342*(7-8), 502-516.

Gischig, V. S., & Preisig, G. (2015, May). Hydro-fracturing versus hydro-shearing: a critical assessment of two distinct reservoir stimulation mechanisms. In *13th ISRM International Congress of Rock Mechanics. OnePetro*.

Gleeson, T., & Ingebritsen, S. (2016). Crustal permeability. (Eds.) Wiley.

Gomila, R., et al. (2021). Quantitative anisotropies of palaeogeomeability in a strike-slip fault damage zone: Insights from micro-CT analysis and numerical simulations. *Tectonophysics, 810*, 228873.

Grawinkel, A., & Stöckhert, B. (1997). Hydrostatic pore fluid pressure to 9 km depth-Fluid inclusion evidence from the KTB deep drill hole. *Geophysical Research Letters*, 24(24), 3273-3276.

Guillou-Frottier, et al. (2020). On the morphology and amplitude of 2D and 3D thermal anomalies induced by buoyancy-driven flow within and around fault zones. *Solid Earth, 11*(4), 1571-1595.

Heap, M. J., et al. (2020). Towards more realistic values of elastic moduli for volcano modelling. *Journal of Volcanology and Geothermal Research, 390*, 10668.

Horne, R. N. (1979). Three-dimensional natural convection in a confined porous medium heated from below. *Journal of Fluid Mechanics, 92*(4), 751-766.

Horton, C. W., & Rogers, F. T., Jr. (1945). Convection currents in a porous medium. *Journal of Applied Physics, 16*(6), 367-370.

Ito, T., & Zoback, M. D. (2000). Fracture permeability and in situ stress to 7 km depth in the KTB scientific drillhole. *Geophysical Research Letters, 27*(7), 1045-1048.

Jiang, C., Wang, X., Sun, Z., & Lei, Q. (2019). The role of in situ stress in organizing flow pathways in natural fracture networks at the percolation threshold. *Geofluids*. https://doi.org/10.1155/2019/3138972.

Jolie, E., Scott, S., Faulds, J., Chambefort, I., Axelsson, G., Guíteŕrez-Negrín, L. C., Regenspurg, S., Ziegler, M., Ayling, B., Richter, A., & Zemedkun, M. T. (2021). Geological controls on geothermal resources for power generation. *Nature Reviews Earth & Environment, 2*(5), 324-339.

Karaoğlu, Ö., Bayer, Ö., Turgay, M. B., & Browning, J. (2020). Thermomechanical interactions between crustal magma chambers in complex tectonic environments: Insights from Eastern Turkey. *Tectonophysics, 793*, 228607.

Karaoğlu, Ö., Bazargan, M., Baba, A., & Browning, J. (2019). Thermal fluid circulation around the Karliova triple junction: Geochemical features and volcano-tectonic implications (Eastern Turkey). *Geothermics, 81*, 168-184.

Karaoğlu, Ö., Browning, J., Bazargan, M., & Gudmundsson, A. (2016). Numerical modelling of triple-junction tectonics at Karliova, Eastern Turkey, with implications for regional magma transport. *Earth and Planetary Science Letters, 452*, 157-170.

Katto, Y., & Masuoka, T. (1967). Criterion for the onset of convective flow in a fluid in a porous medium. *International Journal of Heat and Mass Transfer, 10*(3), 297-309.

Kozlovsky, Y. (1984). The World's Deepest Well. *Scientific American, 251*(6), 98-105.

Lamur, A., et al. (2017). The permeability of fractured rocks in pressurised volcanic and geothermal systems. *Scientific reports, 7*(1), 1-9.

Ledingham, P., Cotton, L., & Law, R. (2019). The united downs deep geothermal power project. In *Proc. 44th Work. Geotherm. Reserv. Eng*, (pp. 1-11).

Liang, X., Xu, T., Feng, B., & Jiang, Z. (2018). Optimization of heat extraction strategies in fault-controlled hydro-geothermal reservoirs. *Energy, 164*, 853-870.

López, D. L., & Smith, L. (1995). Fluid flow in fault zones: Analysis of the interplay of convective circulation and topographically driven groundwater flow. *Water Resources Research, 31*(6), 1489-1503.

Magri, F., et al. (2017). Thermal convection of viscous fluids in a faulted system: 3D benchmark for numerical codes. *Energy Procedia, 125*, 310-317.
Magri, F., Möller, S., Inbar, N., Möller, P., Raggad, M., Rödiger, T., Rosenthal, E., & Siebert, (2016). 2D and 3D coexisting modes of thermal convection in fractured hydrothermal systems-Implications for transboundary flow in the Lower Yar-mouk Gorge. *Marine and Petroleum Geology*, 78, 750–758.

McLellan, J. G., Oliver, N. H. S., & Schaub, P. M. (2004). Fluid flow in extensional environments; numerical modelling with an application to Hamersley iron ores. *Journal of Structural Geology*, 26(6–7), 1157–2117.

Mitchell, T. M., & Faulkner, D. R. (2009). The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: A field study from the Atacama fault system, northern Chile. *Journal of Structural Geology*, 31(8), 802–816.

Morris, A., Ferrill, D. A., & Henderson, D. B. (1996). Slip-tendency analysis and fault reactivation. *Geology*, 24(3), 275–278.

Nur, A. M. O. S., & Walder, (1990). Time-dependent hydraulics of the Earth's crust. In *The Role of Fluids in Crustal Processes Washington, DC* (Vol. 113). National Academy Press.

Ojala, A. E., Mattila, J., Markovaara-Koivisto, M., Ruskeeniemi, T., Palmu, J. P., & Sutinen, R. (2019). Distribution and morphology of landslides in northern Finland: An analysis of postglacial seismic activity. *Geomorphology*, 326, 190–201.

Ord, A., & Oliver, N. H. S. (1997). Mechanical controls on fluid flow during regional metamorphism: Some numerical models. *Journal of Metamorphic Geology*, 15(3), 345–359.

O'Sullivan, M. J., Pruess, K., & Lippmann, M. J. (2001). State of the art of geothermal reservoir simulation. *Geothermics*, 30(4), 395–429.

Parisio, F., Vilarrasa, V., Wang, W., Kolditz, O., & Nagel, T. (2019). The risks of long-term re-injection in supercritical geothermal systems. *Nature communications*, 10(1), 1–11.

Paulillo, A., Cotton, L., Law, R., Striolo, A., & Lettieri, P. (2020). Geothermal energy in the UK: The life-cycle environmental impacts of electricity production from the United downs deep geothermal power project. *Journal of Cleaner Production*, 249, 119410.

Person, M., Hofstra, A., Sweetkind, D., Stone, W., Cohen, D., Gable, C. W., & Banerjee, A. (2012). Analytical and numerical models of hydrothermal fluid flow at fault intersections. *Geofluids*, 12(4), 312–326.

Rowland, J. V., & Sibson, R. H. (2004). Structural controls on hydrothermal flow in a segmented rift system, Taupo Volcanic Zone, New Zealand. *Geofluids*, 4(4), 259–283.

Saevarsdottir, G., Tao, P. C., Stefansson, H., & Harvey, W. (2014). Potential use of geothermal energy sources for the production of lithium-ion batteries. *Renewable Energy*, 61, 17–22.

Sardini, P., Ledésert, B., & Touchard, G. (1997). Quantification of microscopic porous networks by image analysis and measurements of permeability in the Soultz-sous-Forêts granite (Alsace, France). In *Fluid flow and Transport in Rocks*, (pp. 171–189).

Schulz, S. E., & Evans, J. P. (1998). Spatial variability in microscopic deformation and composition of the Punchbowl fault, southern California: Implications for mechanisms, fluid–rock interaction, and fault morphology. *Tectonophysics*, 295(1–2), 223–244.

Scott, S., Driesner, T., & Weis, P. (2017). Boiling and condensation of saline geothermal fluids above magmatic intrusions. *Geophysical Research Letters*, 44(4), 1696–1705.

Sibson, R. H. (1987). Earthquake rupturing as a mineralizing agent in hydrothermal systems. *Geology*, 15(8), 701–704.

Siebenaller, L., et al. (2013). Fluid record of rock exhumation across the brittle–ductile transition during formation of a metamorphic core complex (Naxos Island, Cyclades, Greece). *Journal of Metamorphic Geology*, 31(3), 313–338.

Siler, D. L. H., et al. (2019). Three-dimensional geologic mapping to assess geothermal potential: Examples from Nevada and Oregon. *Geothermal Energy*, 7(1), 1–32.

Sun, Z. X., et al. (2017). Numerical simulation of the heat extraction in EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model. *Energy*, 120, 20–33.

Upton, P., Baxter, K., & O'Brien, G. W. (1998). Coupled mechanical/fluid flow models of trap integrity and fault reactivation: Application to the North West Shelf of Australia. *The APPEA Journal*, 38(1), 488–499.

Violay, M., Heap, M. J., Acosta, M., & Madonna, C. (2017). Porosity evolution at the brittle-ductile transition in the continental crust: Implications for deep hydro-geothermal circulation. *Scientific reports*, 7(1), 1–10.

Wanner, C., Diamond, L. W., & Alt-Epping, P. (2019). Quantification of 3-D thermal anomalies from surface observations of an orogenic geothermal system (Grimsel Pass, Swiss Alps). *Journal of Geophysical Research: Solid Earth*, 124(11), 10839–10854.

Watanabe, N., Abe, H., Okamoto, A., Nakamura, K., & Komai, T. (2021). Formation of amorphous silica nanoparticles and its impact on permeability of fractured granite in superhot geothermal environments. *Scientific reports*, 11(1), 1–11.

Watanabe, N., Egawa, M., Sakaguchi, K., Ishibashi, T., & Tsuchiya, N. (2017). Hydraulic fracturing and permeability enhancement in granite from subcritical/brittle to supercritical/ductile conditions. *Geophysical Research Letters*, 44(11), 5468–5475.

Zareidarmiyan, A., Parisio, F., Makhnenko, R. Y., Salarirad, H., & Vilarrasa, V. (2020). How equivalent are equivalent porous media? *Geophysical Research Letters*, 48, e2020GL0089163.

Zoback, et al. (2003). Determination of stress orientation and magnitude in deep wells. *International Journal of Rock Mechanics and Mining Sciences*, 40(7–8), 1049–1076.

Zoback, M. L. (1992). Stress field constraints on intraplate seismicity in eastern North America. *Journal of Geophysical Research: Solid Earth*, 97(B8), 11761–11782.