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

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.


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 meterwide geological structures that localize deformation (Ben-Zion & Rovelli, 2014) and modify the petrophysical properties of the crust from the surface down to the brittle-ductile transition (BDT). CFZs are defined by faults intersection and interconnected 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., 2017Watanabe et al., , 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 heatproducing 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;. 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 affected by stresses, whether on their pressure (undrained 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 vertical 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;. Thus, the stress orientation has an effect on fluid flow, as already 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 aspects 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 regimes 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 faulted 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 another and will impact the formation of a high-temperature 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 wellknown tool offering a complete access to the solution of partial differential equations. Benchmark tests have already been performed with previous numerical codes (OpenGeoSys, Comsol Multiphysicsä 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 temperature of 20°C and a pressure of 10 5 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 À2 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 assumed to be pure water, and the fluid density depended 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, hydraulic and mechanical parameters are detailed in Table 1.
The Young's moduli imposed in numerical models 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.
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 v ), maximum horizontal (S Hmax ) and minimum horizontal (S hmin ) components, which are used regularly in geomechanical studies of reservoirs (Anderson, 1905;Zoback et al., 2003): The relative 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.
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: 3 = 1 × 10 6 + (15000 × ( )) 2 = 4 × 10 6 + (20000 × ( )) 1 = 7.5 × 10 6 + (28667 × ( )) ð1Þ Low stress intensity: where r 1 , r 2 ; r 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., r 1 ¼ S Hmax , r 2 ¼ S v , r 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: where x represents the degree of rock fracturing; x ¼ 1 for fully fractured rock and x ¼ 0 for fully intact rock. Here, we took x ¼ 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 for fractured rock. The k i is permeability (m 2 ) of intact rock, and k f is permeability (m 2 ) 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: where s is shear stress and r 0 n is effective normal stress defined as: where r n is total normal stress; a B is the Biot-Willis coefficient and pf is fluid pressure.

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 9 10 -9 m.s À1 at the center of the fault and 20 9 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  (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 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. 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 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.
( Figs. 2 and 3). The amplitudes of the second temperature anomaly were + 35°C and + 40°C for low and high stress intensity, respectively. In the horizontal cross section, the two temperature anomalies were found for each intensity. In high stress intensity, 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 limited. 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 temperature 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 anomalies spread largely beyond the fault.
Temperature anomalies of + 25°C were found in the basement, suggesting that in a strike-slip system, the positive temperature anomaly, and independently 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. Finally, in compressional regime these anomalies were localized in the near vicinity of the fault. The tectonic regimes have shown to play a key role in temperature distribution, and this was clearly related 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 movements at the ends of the fault (Figs. 2 and 3). In high stress intensity the maximum fluid velocity was 30 9 10 -9 m.s À1 , against 21 9 10 -9 m.s À1 in low stress intensity. The minimum fluid velocities were, for each intensity, 1 9 10 -9 m.s À1 . 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 9 10 -9 m.s À1 , compared to 17 9 10 -9 m.s À1 in low stress intensity. The minimum fluid velocity was 5 9 10 -9 m.s À1 in high stress intensity, compared to 1 9 10 -9 m.s À1 in low stress intensity.
In compressional tectonic regimes, for each stress intensity, there was a slightly different convective pattern. In high stress intensity, there were two upward and two downward movements, whereas in low stress intensity there were two upward and three downward movements. The maximum fluid velocity for high stress intensity was 5 9 10 -9 m.s À1 and for low stress intensity was 10 9 10 -9 m.s À1 . For this tectonic regime, the maximum fluid velocity was present when the lowest stress intensity was applied. In horizontal cross section, at 2 km depth and for high stress intensity, upward movements were localized where temperature anomalies were positive, whereas downward movements were localized where temperature anomalies were negative. The fluid velocity was between 1 and 2 9 10 -9 m.s À1 . In low stress intensity, the downward movements were localized at the extremities and at the center of the fault. The upward movements were localized between each downward movement. The maximum fluid velocity was 10 9 10 -9 m.s À1 and the minimum was 1 9 10 -9 m.s À1 . At a depth of 2 km, fluid velocities varied from 1 to 2 9 10 -9 m.s À1 . 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 9 10 -9 m.s À1 and the maximum was 12 9 10 -9 m.s À1 at the downward movements. In high stress intensity, the minimum fluid velocity was 1 9 10 -9 m.s À1 and the maximum was 17 9 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. 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 steadystate regime at 203 kyrs for high stress intensity and 223 kyrs for low stress intensity.
The benchmark experiment reached the steadystate 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 9 10 -12.1 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 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 Figure 6. Slip tendency in high stress intensity application (in red), and low stress intensity application (in green). 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 strikeslip 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 (Karaog 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, that localize fluids flow (Person et al. 2012;Karaog lu et al., 2016Karaog lu et al., , 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.

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 Karaog 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://creativecom mons.org/licenses/by/4.0/.

APPENDIX 1
The calibration of our numerical experiments was performed based on results from Comsol Multiphysicsä version 3.5  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 9 5.5 9 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 9 10 -15 m 2 . The fluid flow velocity was 2.03 9 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 9 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 Dar-cyÕ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).