Advanced 3D TH and THM Modeling to Shed Light on Thermal Convection in Fault Zones With Varying Thicknesses

Fault zones exhibit 3D variable thickness, a feature that remains inadequately explored, particularly with regard to the impact on fluid flow. Upon analyzing an analytic solution, we examine 3D thermal‐hydraulic (TH) dynamical models through a benchmark experiment, which incorporates a fault zone with thickness variations corresponding to realistic orders of magnitude. The findings emphasize an area of interest where vigorous convection drives fluid flow, resulting in a temperature increase to 150°C at a shallow depth of 2.7 km in the thickest sections of the fault zone. Moreover, by considering various tectonic regimes (compressional, extensional, and strike‐slip) within 3D thermal‐hydraulic‐mechanical (THM) models and comparing them to the benchmark experiment, we observe variations in fluid pressure induced by poroelastic forces acting on fluid flow within the area of interest. These tectonic‐induced pressure changes influence the thermal distribution of the region and the intensity of temperature anomalies. Outcomes of this study emphasize the impact of poroelasticity‐driven forces on transfer processes and highlight the importance of addressing fault geometry as a crucial parameter in future investigations of fluid flow in fractured systems. Such research has relevant applications in geothermal energy, CO2 storage, and mineral deposits.


Introduction
In light of global change, understanding the dynamics of current and/or past fluid flows within faults in the solid Earth is of paramount scientific and economic significance.On a worldwide scale, fault zone are ubiquitous in various tectonic regimes (Scibek, 2020).Fault zones are intricate volumes of fractured, altered and deformed rocks with variable geometry, coevolving across three spatial dimensions and one temporal dimension (Wibberley et al., 2008) and playing a central role in various processes that have both fundamental (e.g., deformation localization) and far-reaching societal impact such as geothermal energy, CO 2 storage, and mineral deposits (Handy et al., 2007;Jolie et al., 2021;Nicol et al., 2017;Phillips, 1972).Despite numerous studies related to these applications, characterizing and predicting fluid flows within these systems remains a challenge.Fluid flow depends on the interaction between Thermal (T), Hydraulic (H), Mechanical (M) and Chemical (C) parameters (Ingebritsen & Appold, 2012;Rutqvist et al., 2002).Using TH and THM coupling, this study aims to investigate the impact of a Geometrical (G) fault parameter on the dynamic behavior of fluid flow.
Fault zones are complex structures characterized by a succession of undeformed protolith, damage zones, and fault cores that reflect intense deformation (e.g., Faulkner et al., 2003).These damage zones and fault cores contain gouges, breccias, cataclasite relay branch lines, bends, and may exhibit unusually high permeability values (Bischoff et al., 2024;Gleeson & Ingebritsen, 2016).Deep and anomalously permeable zones, known as Crustal Fault Zones (CFZ), are influenced by the mechanical heterogeneity of the surrounding rock (Childs et al., 1996;Ferrill & Morris, 2003;Peacock & Sanderson, 1992).The irregular and segmented propagation of these faults within a rock volume depends on factors such as the amount of displacement they accommodate (Childs et al., 1995(Childs et al., , 2009;;Cox & Scholz, 1988;Ferrill & Morris, 2001;Ferrill et al., 1999;Manighetti et al., 2004;Naylor et al., 1986;Perrin et al., 2016).For fault zones, the direction, length, depth, dip, and thickness constitute the set of parameters defining the geometrical (G) position of the fault in space.In 1875, in Colorado (USA), J. W. Powell observed different thicknesses within a single fault zone, including a thin deformation zone (several centimeters) without notable filling and a thicker deformation (several tens meters) zone characterized by block filling (Powell, 1875).On another illustrative instance, in the Upper Rhine Graben (URG) of Western Europe, fault zone thickness is closely related to a pre-structuring of Variscan tectonics, which has been reactivated during polyphase deformation associated with Mesozoic and Cenozoic graben tectonics (Valley, 2007).Several geothermal boreholes have intersected naturally permeable faults developed in Carboniferous monzogranite rocks (Vidal & Genter, 2018).In the French part of the URG, these faults are exploited by the geothermal industry for electricity production at Soultz-sous-Forêts and high-temperature heat purposes at Rittershoffen (Dalmais et al., 2022).A detailed structural analysis of one exploration borehole revealed a cluster organization of the 3,400 fractures observed along 1,200 m of core length (Genter & Traineau, 1996).Fracture density within the Carboniferous monzogranite is 1.5 fractures/cm and fault thicknesses derived from core observations in the monzogranite range from centimeters to multiple decimeters (Genter & Traineau, 1996).In the basement, the highest concentrations of fractures are associated with fault zones, which are closely related to hydrothermal alteration due to past and present fluid circulations.When different faults intersect, fault zones can easily reach thicknesses of 10 to hundreds of meters and even up to several kilometers (Ben-Zion & Rovelli, 2014).For instance, in Andalusia (Spain), the Carboneras vertical strike-slip fault zone is described with a thickness of 1 km (Keller et al., 1995;Rutter et al., 1986).However, where the Carboneras Fault connects to the Palomeras Fault, the described fault zone has a thickness of about 4 km (Faulkner et al., 2003).Seebeck et al. (2014) demonstrated that fault zone thickness can vary by three orders of magnitude over a "short distance", potentially directly impacting fluid flow.
Fault zone might be underlined by past and/or present fluid flow.Past fluid flow within branch lines and relay zones has been inferred from Zn-Pb deposits in the Irish orefield (Carboni et al., 2003) and the Pontgibaud Crustal Fault Zone (Bouladon et al., 1961). Furthermore, Garden et al. (2001) demonstrated that these CFZ serve as migration channels for hydrocarbon systems.In these naturally permeable and deep systems, present-day fluid flow is acknowledged (McCaig, 1988), and convective heat transfer has been suggested as a potential driving force (Etheridge et al., 1984;Fyfe, 2012;Horne, 1979).Without an additional heat source (e.g., a magma chamber), the presence of fluid in an anomalously permeable and deep environment allows for the emergence of a positive thermal anomaly.This trend has been noted in previous studies (Duwiquet, 2022;Moeck, 2014), but without considering variations in fault zone thickness, which are present in some cases (e.g., Carboneras Fault Zone in Andalusia, Spain, Faulkner et al., 2003); Margeride Fault in the French Massif Central, Talbot et al., 2005).Malkovsky and Pek (2004) demonstrated that, among other factors, fault zone thickness is a key control parameter for heat transfer modes.Malkovsky and Magri (2016) examined these effects on flow regimes within a vertical fault zone, by considering the viscosity variation with temperature.Futhermore most of numerical studies of fluid flow in faulted zones do not consider a variable fault thickness and the impact of tectonic regimes has not been specifically tackled (see however Duwiquet et al., 2021b;Duwiquet, Magri, et al., 2022;Guillou-Frottier et al., 2024).Here, we propose to investigate numerically these effects in a fault zone with variable thickness, and with a realistic fluid density law.
Various forces within fault zones impact fluid flow.Buoyancy-driven force can result in thermal convection due to temperature-dependent fluid density.Other factors, such as topography (Forster & Smith, 1989) and metamorphic or magmatic fluid production (Etheridge et al., 1984), create pressure gradients that affect fluid flow.The poroelasticity-driven force, describes the interaction between fluids and deformation in the porous medium.It has also been identified as a source of pressure variations that directly impacting fluid circulation.Duwiquet, Magri, et al. (2022) and Duwiquet et al. (2024b) have studied the poroelasticity-driven force in a faulted system with a thickness of 0.4 km, revealing different convective dynamics for each tested tectonic regime (extensional, compressional, strike-slip).However, these effects have not been examined in larger fault zones, and even less in fault zones with variable thickness.

Materials and Methods
The Comsol Multiphysics™ software, based on the Finite Element Method, can couple various time-dependent physical processes, such as fluid flow, heat transfer, and elastic deformation of materials, within a 3D geometry.The fault zone has a thickness ranging from 0.4 to 2 km, with a depth of 5 km (Figure 1).The shape of the fault zone is considered as a simplified shape, characterized by a parallelepiped, in such a way as to be able to evaluate in 3 dimensions the impact of a thickness variation on fluid flows.Within the fault zone, permeability varies spatially (Appendix A).This variation corresponds to lithologic changes between damaged and core zones along the width and length of the fault zone.This permeability variation, corresponding to the conceptual "multiple fault core" model, has already been demonstrated on the Carboneras fault (Faulkner et al., 2003(Faulkner et al., , 2010)).This conceptual model can be applied both to large fault zones (Duwiquet et al., 2019(Duwiquet et al., , 2021b) ) and also to models simplified fault zones (Duwiquet, Guillou-Frottier, et al., 2022).Furthermore, considering the compaction effects, permeability decreases with depth according to Ingebritsen and Manning (2010).By integrating both a permeability decrease with depth (z) and lateral variations (x, y), the heterogeneous permeability can be written as: As for all tectonic regimes analyzed (compressional, extensional, strike-slip tectonic regimes), the stress intensity increases positively with depth.For the compressional system, S Hmax is applied perpendicular to the fault zone, S hmin is applied on the faces parallel to the fault zone, and S v is applied vertically.
where K f is the fixed space-dependent permeability of the fault and K f0 × 201 [m 2 ] represents the maximum permeability value at the surface.The length δ (m) characterizes the intensity of the decrease in permeability with depth.Considering the system as intensely fractured, and in accordance with the models of Garibaldi et al. (2010), Guillou-Frottier et al. (2013) and Duwiquet et al. (2019) the chosen value for α is 2,500 m.Opting for constants 100 and 101 allows for the fixed variation of permeability across two orders of magnitude.The last term on the right side of Equation 3 causes the permeability to alternate along a sinusoid applied to the x, y, and z axes.The term λ corresponds to the wavelength of the sinusoid.To reproduce relatively fine alternations of high and low permeability, as suggested by field observations, a low value of λ is chosen.
The fault zone is located in the middle of a 14 × 14 × 7 km basement (Figure 1).These dimensions have been pretested to avoid any potential edge effects occurring within and around the fault zone.The numerical modeling results are presented in steady-state.A no-flow condition is imposed on all boundaries except on the upper one, for which an atmospheric pressure of 10 5 Pa and a temperature of 10°C are fixed (Figure 1).The external vertical faces are adiabatic.The initial thermal regime corresponded to a geothermal gradient of 30°C/km, and equivalent to the average geothermal gradient of the European continental crust.Then, a heat flux of 100 mW/m 2 is imposed at the lower boundary (Figure 1).This heat flux represents the combined effects of heat flux from the mantle and the heat emitted by the decay of radioactive elements in the crust, and is consistent with measured surface heat flow in the French Massif Central (Lucazeau et al., 1984).These boundary conditions are imposed on the whole domain (both the fault zone and the basement).Based on these boundary conditions, the evolution of a conductive heat transfer regime toward a convective regime, where the fluid flows are seen along the streamlines, is examined (Figure 2).The fluid is composed of pure water.The fluid density depended on pressure and temperature conditions (as detailed in Duwiquet et al., 2021b).The couplings presented here are TH couplings for the benchmark experiment (Figure 2) and THM couplings for the tectonic regimes tested (Figures 5-7).THM couplings operate through pressure and temperature variables.Temperature is solved by the heat equation but also appears in the density and viscosity laws, both parameters being present in the Darcy law.For the mechanical coupling, based on the poroelasticity assumption, the fluids in a reservoir are affected by stresses, either on their pressure (undrained conditions in low-permeability media with, for example, an increase in pressure in a state of compressive stress), or on their circulation (drained conditions in permeable media with, for example, convection from the most compressed to the least compressed regions).Fluid flow is examined for each tectonic regime in the THM models.The consequences in terms of temperature anomalies ΔT [°C] are also described.The temperature anomalies ΔT, correspond to abnormally hot and abnormally cold temperatures compared to the undisturbed initial conductive thermal regime with the 30°C/km temperature gradient.
Details of the coupling between Darcy's law, the heat equation, and Hooke's law can be found in Duwiquet et al. (2021b).The mesh is defined by 57,841 tetrahedra, with maximum mesh sizes of 400 m for low permeability and 80 m for high permeability.A preliminary convergence test showed that a finer mesh gave the same results.Futhermore, the vertical boundaries effects are limited as long as we look at the fluid flows where the permeability is the higher, that is, in the fault zone.All parameters integrated into this model are provided in Table 1.Each parameter is defined as a function of the domain (i.e., fault zone or basement).The transition between a parameter fixed in the fault zone and the same parameter in the basement is sharp.
In order to assess the impact of tectonic regimes on fluid flow, we propose comparing the results of numerical experiments with a benchmark experiment.Here, the benchmark experiment considers a TH coupling.By maintaining the same geometry and physical properties, different tectonic regimes are implemented.For these models, we have free mechanical boundary conditions 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 are applied.We aim to study the impact of the three main tectonic regimes, namely compressional, extensional, and strikeslip.We thus create three numerical models, one for each tectonic regime.For this, we consider an Andersonian assumption-where the principal stresses are expressed with vertical S V , maximum horizontal S Hmax , and minimum horizontal S hmin components -which is regularly used in geomechanical studies of reservoirs (Anderson, 1905;Zoback et al., 2003).The relative magnitudes of these stresses determine the simplified tectonic regime model: Assuming that our models, align with the principal stresses, pure normal stresses are applied on the lateral boundaries (no shear applied on the boundaries).For the compressional regime, the fault is perpendicular to the maximum horizontal stress S Hmax applied on the boundaries (Figure 1).For the extensional regime, the fault is perpendicular to S hmin .Lastly, for the strike-slip regime, the fault is at 45°between S hmin and S Hmax orientation.Vertical stress S V can be expressed as a function of the overlying rock mass and then adjust the horizontal/vertical stress ratios to accommodate different tectonic regimes: where ρ s [kg/m 3 ] is the rock density (see Table 1), g[m/s 2 ] is the acceleration of gravity, and z[m] is the vertical upwards axis.An increase in depth (z < 0) results in more compressive stress (positive compression convention).α Hmax [ ] and α hmin [ ] are the horizontal-to-vertical stress ratios.To accommodate different tectonic regimes, we must set: The values for the stress ratios remain to be estimated.A comprehensive parametric study would require testing various values for each tectonic regime.However, to follow a step-by-step approach and limit the number of unknowns, we decided to fix these values based on literature data.The San Andreas stress regime, undergoing a strike-slip regime, is selected as it is a well-documented system and results in a highly deviatoric configuration, thus intensifying the stress effects (Zoback, 1992).
where stresses are expressed in Pa.
Arguing that the static parts become less significant at deep depths (several kilometers), we focus on the gradients in Equation 3to adapt the stress ratios to our models.Based on the gradient values in Equation 3 and the principal/ Andersonian relations in Table 2, we can propose the set of values presented in Table 3.
It is now possible to investigate the effect of tectonic regimes while maintaining the realistic stress ratios given in Table 3 and illustrated in Figures 1b  and 1c.To facilitate the computation of early solutions with this THM coupling, the stress application is progressively applied from t = 0 years until t = 10 years.

Results
In order to observe any potential impact of the poroelasticity-driven force and the fault thickness as a geometrical (G) parameter, we will first detail the benchmark experiment results.The characterization of fluid flow and thermal distribution within the faulted system will be detailed, alongside an assessment of fluid pressure variations and an analysis of the critical Rayleigh number.Following this, we will describe the results of numerical modeling considering the different tectonic regimes tested in the same manner.A comparison between these results will then be made, allowing us to discuss the impacts of both geometry and poroelasticity-driven forces on fluid flow within naturally fractured systems.

Results of 3D Modeling and Analytical Solution: Benchmark Experiment (TH Coupling) and Critical Rayleigh Number Analysis
The initial results of 3D numerical modeling (TH coupling) represent the benchmark experiment of our study.Without stress application, fluid flow takes place on both sides of the fault (Figure 2).As the fluid density increases, its sinks within the fault volume, leading to a convection pattern related to buoyancy forces.Under these conditions, this free convection generates a thermal disturbance.In initial regime (i.e., in purely diffusive setting), the 150°C isotherm is set at a depth of 5 km, whereas, in the area with the thickest deformation zone, the same isotherm is found at a depth of 2.7 km.In the thinnest zone, the 150°C isotherm rises to a depth of 3.3 km (Figure 2).In the thickest part of the fault, fluid flow is defined by a vigorous convective cell.The wavelength of the convective cell occupies only one-third of the fault surface area (Figure 2).Fluid velocities vary between 1 × 10⁻⁹ m/s and 16 × 10⁻⁹ m/s.The lowest velocities are concentrated in the bottom part of the convective cell, while the highest velocities are concentrated in the top part of the convective cell.In this benchmark experiment, the convective cell is located in the center of the thickest part of the fault zone, and it exhibits a typical cellular pattern.The results will differ after applying different tectonic regimes.
In the thinnest part of the fault, fluid flow is also apparent (Figure 2).Less powerful, the fluid upwelling is barely visible at the bottom of the fault.The upwelling of the isotherms is more indicative of this convection, which, while less strong, has a significant wavelength, as it occupies two-thirds of the fault surface area.
Clearly, the thickness of the fault zone affects the vigor of convection.It seems that the onset of convection would occur for a critical thickness, meaning that the critical Rayleigh number (above which thermal convention occurs) depends on fault thickness.Indeed, Malkovsky and Magri (2016) studied the onset of thermal convection within a permeable faulted zone of finite width through a linear stability analysis.They improved upon the preliminary study by Malkovsky and Pek (1997), where viscosity was assumed to be constant.The fault zone has a halfthickness δ, a height H, with a fluid viscosity that is temperature-dependent.They showed that the critical Rayleigh number is expressed as: where a is a constant equal to 1 for the constant viscosity case, and 2.466 for a temperature-dependent viscosity with an global average temperature gradient (see Malkovsky & Magri, 2016).Figure 3 illustrates the role of the fault width on the critical Rayleigh number for different fault heights.For wide fault zones (hundreds to thousands of meters), the critical Rayleigh number is decreased by 1-2 orders of magnitude (compared to a fault width of tens of meters), and thermal convection occurs more easily.
For a small fault thickness, the upwelling of fluids is less significant, and the wavelength is greater.Conversely, for a large thickness, the convection is more intense, and the wavelength is smaller.This confirms the result of the benchmark experiment (Figure 2).In order to see if the tectonic regimes can impact the fluid flow and to assess the potential effect on the thermal distribution of the system, we will focus on the area where the fluid flow by convection is effective in steady-state, in the thickest zone (Figures 5-7).
Then, we will investigate the fluid pressure variation in this area.
In the thickest part of the fault zone in the benchmark experiment, the fluid pressure remains constant at around 19.6 MPa at a depth of two km (Figure 4).Without stress application, the pressure here corresponds to the hydrostatic pressure.A slight difference within the fault can be explained by a fluid passage with high velocity (15 × 10 9 m/s).In this case, with no other  significant variations in the pressure gradient known, such as those related to topography (which is flat in this case) and/or the production of metamorphic or magmatic fluids, the convection is considered to be free.Taking a poroelastic hypothesis into account and after considering different tectonic regimes, the fluid pressure will vary, causing changes in the general convective dynamics.

Compressional Tectonic Regime
In the compressional tectonic regime (Figure 5), the fluid flow pattern is characterized by downward and upward movements along the Z and X axes.Fluid flow velocities range from 1 × 10 9 m/s to 10 × 10 9 m/s (Figure 5).The fastest fluid flows are located in the top part of the upward movement, while the slowest fluid flows are located in both the downward movement and at the beginning of the upward movement.Where the fluid motion changes from an upward to a downward movement at the top of the convective cell, the cell is "pinched," causing the typical cellular pattern (Figure 2) to become distorted.In this scenario, the fluids almost go from the highest to the lowest velocity recorded, from 9 × 10 9 m/s to 1 × 10 9 m/s, within a small volume.The overall fluid velocities are less significant, and the pinching at the top of the convective cell has not been previously observed.The upward fluid flow is accompanied by the red zone, which represents the positive temperature anomaly, with a maximum value of 47°C (Figure 5).At the bottom of the convective cell, the pinched shape is also found in the positive temperature anomaly.Here, unlike in the benchmark experiment, the maximum temperature anomaly is found where the fluid flow velocity is the fastest.Thus, in comparison, there are differences in terms of fluid flow velocities, thermal distributions within the system, and the presence of a pinched shape at the top of the convective cell.
In the compressional tectonic regime, fluid pressure in both the basement and the fault zone is higher than in the benchmark experiment (Figure 5).The fluid pressures in the basement and fault zone are 23 and 21.25 MPa, respectively (Figure 5).This difference of 1.75 MPa leads the fluids toward the zone of lower pressure and higher permeability, that is, the fault zone.This difference may impacts the general convective dynamics in the fault zone, where a distortion of the convective cell is observed and has consequences on the thermal distribution of the system (Figure 5).Since differences exist between the benchmark experiment and the compressional tectonic system, the buoyancy-driven force alone cannot explain the current convective patterns.After stress application, the poroelasticity-driven force, coupled with the buoyancy force, both have an impact on fluid circulation.

Extensional Tectonic Regime
In the extensional tectonic regime (Figure 6), the fluid flow pattern is characterized by both downward and upward movements, primarily along the Z and X axes (Figure 6).Notably, fluid movements along the Z, X, and Y axes are present in the upper part of the cell.Fluid flow velocities range from 0.8 × 10 9 m/s to 17.7 × 10 9 m/s.These velocities are of the same order of magnitude as the fluid velocities in the benchmark experiment and are faster than those in the compressional tectonic regime.However, compared to the compressional tectonic regime, the fastest velocities remain in the upward movement of the fluid, and the slowest velocities are in the downward movement.When the fluid motion switches from an upward to a downward movement at the top of the convective cell, the "pinching" previously described is more pronounced.Overall, the convective cell is more distorted than in the previous cases, resulting in a fluid flow initiated in all three spatial dimensions (Z, X, Y axes).
The positive temperature anomalies follow the upward and fastest fluid movements (Figure 6).The maximum positive temperature anomaly, 43°C, is at the top of the convective cell, where the convective cell is pinched.The general shape of the positive temperature anomaly follows the fluid flow distortion.In this case, the pinching of the convective cell, more pronounced than in the compressional tectonic regime, facilitates the development of the positive temperature anomaly along the Z and Y planes.
In the extensional tectonic regime, the fluid pressure, both in the basement and the fault zone, is lower than that in the benchmark experiment (Figure 6).The fluid pressure difference between the basement and the fault zone is 1.2 MPa, with a fluid pressure in the basement of 19.5 MPa, which decreases until it reaches 18.3 MPa in the fault zone.Similarly, as in the compressional tectonic system, the difference in fluid pressure between the basement and the fault leads the fluid toward the fault zone, defined by lower pressure and higher permeability.This difference may impacts the general convective dynamics in the fault zone, where a distortion of the convective cell is observed and has consequences on the thermal distribution of the system (Figure 6).As differences exist between the benchmark experiment and the compressional tectonic system, the buoyancy-driven force alone cannot explain the convective cells present in the thickest part of the fault zone.After stress application, the poroelasticity-driven force, coupled with the buoyancy force, both impact fluid circulation, as in the compressional tectonic regime.

Strike-Slip Tectonic Regime
In a strike-slip tectonic regime (Figure 7), fluid flow is characterized by continuous upward and downward movement along all three axes (Z, X, Y).Fluid flow velocities range from 1 × 10 9 m/s to 13.8 × 10 9 m/s.The fastest velocities are concentrated in the upward movement, while the slowest velocities are focused on the downward movement.These velocities are within the same orders of magnitude as in the previous experiments.The pinching, previously described in compressional and extensional tectonic regimes, gives way to a continuous fluid motion developed across all three spatial dimensions.This convective dynamics has already been observed in other fault/basement systems modeled (Duwiquet et al., 2021a(Duwiquet et al., , 2021b) ) and has been identified as a "doublelike convective pattern."As described earlier, the fastest upward movement of fluids follows the positive temperature anomaly.The maximum temperature anomaly obtained is +73°C (Figure 7), which is the most intense positive temperature anomaly observed among the different tectonic regimes tested.The upper parts of the double-like convective pattern concentrate the warmest zones.The shapes of the anomalously hot zones, which follow fluid flow, propagate in all three spatial dimensions.
The fluid pressure variation follows the same trend as in the two previous tectonic regimes: a decrease in fluid pressure as the fault zone is approached.However, unlike the other two tectonic regimes (i.e., compressional and extensional), the difference in fluid pressure between the fault and the basement is 0.9 MPa (Figure 7).This represents the lowest fluid pressure difference value of the two other tectonic regimes tested.Compared to the benchmark experiment, the fluid pressure is above the reference value of 19.6 MPa.Thus, as in the previous tectonic regime tested, the fluid is led toward the low-pressure, high-permeability fault zone.However, unlike the other tectonic regimes tested, the small difference in fluid pressure between the basement and the fault zone could allow for more pronounced convective dynamics in all three dimensions of space.Since differences exist between the benchmark experiment and the compressional tectonic regime, the buoyancy-driven force alone cannot explain the convective cells present in the thickest part of the fault zone.After stress application, the poroelasticity-driven force, coupled with the buoyancy force, both impact fluid circulation and have consequences on the thermal distribution of the medium.

Discussion
Geosciences are emerging as a key element in the effective transition from fossil fuels to low-carbon energy and strategic minerals associated with it.In the Earth crust, anomalously permeable areas, such as fault zones, are attractive targets for geothermal energy, CO 2 storage and mineral exploration.It is widely recognized that fluid flows in the crust are influenced by full Thermal, Hydraulic, Mechanical and Chemical parameters as well as various driving forces.These preliminary results demonstrate that the variation of fault thickness as a Geometrical (G) parameter and the poroelasticity-driven force both impact the convective dynamics of fluid flow within these anomalously deep permeable systems.

Poroelasticity Driven Force, One of the Forces That Drive the Fluid
The forces affecting the fluid are diverse and vary in importance.According to Darcy's law (Darcy, 1857), fluid flow behavior is related to the pore fluid pressure gradient: pore fluid always flows from areas of high to low pressure (Terzaghi, 1925).Our numerical modeling shows that the poroelasticity-driven force, generated by the tectonic regimes, alters the general convective dynamics within the thickest deformation zone.Without any topography and without magmatic or metamorphic fluid production, fluid flow results from both buoyancy and poroelastic-driven forces.
In the benchmark experiment, fluid velocities are faster than in other experiments (Figures 2, 5-7).Moreover, in the benchmark experiment, the highest velocities are found in both upward and downward fluid movement (Figure 2).In contrast, in other tectonic regimes, the fastest fluid velocities are found in the upward movement (Figures 5-7).This is unlike the gravitational force, which has a more significant impact at low velocities.However, considering a constant thickness fault zone geometry at 0.4 km, the results of Duwiquet, Guillou-Frottier, et al. (2022) and Duwiquet, Magri, et al. (2022) show a different trend.Indeed, the fastest fluid velocities are located in downward movements.The impetus generated by the difference in fluid pressure between the fault and the basement, following the application of tectonic stresses, along with buoyancy forces favoring the upwelling of less dense hot fluid, could explain this effect, which is more visible in larger fault zones.
Although the thickness of the fault zone depends on many factors, all authors describe heterogeneity and/or anisotropy within each of these structures.In this study, we consider permeability as heterogeneous within the fault zone in all three spatial dimensions.For fluid flow, it is well accepted that heterogeneity and anisotropy are closely related properties.Inhomogeneous materials generally appear homogeneous but anisotropic when considered on a larger scale than heterogeneity.An important consequence is that anisotropy, like heterogeneity, depends on the study scale (Dagan, 1986).

What Impact on Geothermal Potential and Mineral Resources?
Numerical modeling investigation allows us to test the effect of geological and physical parameter variations and to account for fluid flow on the thermal distribution of the system.Considering a 2D model with TH coupling, Duwiquet et al. (2019) highlighted the interest of targeting subvertical/vertical abnormally deep permeable zones for high-temperature geothermal energy.Although conceptually proven, this trend is found in complex natural systems.For example, since 2020, the United Downs Deep Geothermal project has been attempting to target the Porthtowan subvertical Fault Zone reservoir within the Cornish granite (heat producer) (Ledingham et al., 2019;Paulillo et al., 2020).Another example is the Tocomar geothermal system (Argentina), which is controlled by the Calama-Olacapato-El Toro subvertical fault zone (Filipovich et al., 2022).In this system, the positive temperature anomaly is much higher (Figure 7) for a strike-slip tectonic regime.
Figure 3 shows that, without tapering, convection occurs more easily in fault zones characterized by several hundred to several thousand meters thick than in fault zones of 10 of meters thick.When the thickness is variable, the results emphasize also the interest of targeting the thickest parts of fault zones (Figures 2 and 5-7).In 3D, isotherms ascend concurrently with the upward fluid movement.However, the substantial permeability of the surface induces a downward trend (Guillou-Frottier et al., 2020) and leads the emergence of a negative temperature anomaly.These results could be considered and compared with natural systems.For example, in the Margeride fault zone (French Massif Central), its thickness is 1,500 m at the Mentières level, while a few kilometers to the southeast, this thickness is only 400 m.Additional favorability criteria can be considered in the same way as the thickness of the fault zone, such as the amount of displacement (Guillou-Frottier et al., 2024).
Thermally driven convective cells allow the transport of minerals, such as uranium (Li et al., 2021).The surface thickness of the Wollaston-Mudjatik Transition Fault Zone (Rae-Hearn craton -Western Canadian Shield) varies between 11 and 25 km (Poh et al., 2022).It can be observed that where the thickness is largest (>20 km), there is the greatest number of occurrences and/or deposition of uranium (e.g., McClean Lake, Eagle Point, and Rabbit Lake).Within the URG, it is well-known that natural fluid convection is related to faults deeply rooted within the late Carboniferous crystalline fractured basement.Beyond heat and power production, understanding permeability related to fault thickness has further economic applications, such as lithium extraction from geothermal brines circulating in those fault zones.It has been proven that it is possible to extract geothermal lithium during the exploitation of the Soultz power plant (Fries et al., 2022).As the geothermal fluid is mainly circulating in fault zones, understanding the role of fault thickness becomes a relevant challenge to consider for lithium, uranium, and related exploitation.

Conclusions and Perspectives
Fault thickness, as a Geometrical fault parameter (G), plays a key role in controlling fluid flow in the Earth's crust.After analyzing the critical Rayleigh number and employing a simplified numerical modeling approach, these results demonstrate that thermal convection occurs more efficiently in the largest part of the fault zone, but also that the highest positive thermal anomaly is located at 1.5-2.0km from the largest zone.The poroelasticity-driven force, inferred from the tectonic regimes, modifies the general convective dynamics within this area of interest, impacting both the thermal distribution of the system and the potential associated mineralization.At a time of global change, this methodology and the associated results highlight some key parameters to consider during the exploratory phase for heat and/or electricity production, as well as any related mineralization within areas of fundamental scientific and economic interest, such as fault zones.3) corresponds to lithological variations (i.e., fault core/damage zone/basement) in the "multi-core fault" conceptual model (Faulkner et al., 2003(Faulkner et al., , 2010)).The chosen permeability values align with laboratory measurements conducted according to the Heap and Kennedy protocol (2016) on samples collected from the Pontgibaud crustal fault zone in the French Massif Central (Duwiquet et al., 2021b).Furthermore, these values are in line with the permeability data found in the Scibek (2020) database.

Figure 1 .
Figure 1.Model setup and boundary conditions: Meshes are defined by tetrahedrons.The mesh is refined for the highest permeable zone, in this case, the fault (a).The contrast between a refined mesh in the fault and a fine one in the basement causes a numerically coarse meshing at the fault edges.(b) and (c) are stress application examples for a compressional tectonic regime.The length of the red arrows is proportional to the applied stress intensity.As for all tectonic regimes analyzed (compressional, extensional, strike-slip tectonic regimes), the stress intensity increases positively with depth.For the compressional system, S Hmax is applied perpendicular to the fault zone, S hmin is applied on the faces parallel to the fault zone, and S v is applied vertically.

Figure 2 .
Figure 2. 3D model results of the benchmark experiment (i.e., without tectonic regime implementation) are given in steady-state.The highest temperatures are redcolored, and the lowest temperatures are blue-colored.In the initial conductive regime, the 150°C isotherm is located at a depth of 4.6 km.In steady-state regime, fluid flows is marked by solid lines, with direction indicated by colored arrows.Fluid velocities are shown, with red representing the fastest and blue the slowest.Powerful convection is most pronounced in the thickest zone, where the 150°C isotherm (white dotted line) rises to a depth of 2.7 km (see text for mor details).Figures 5-7, with stress application, will focus where fluid flows is effective in the steady-state regime: in the largest part of the fault.

Figure 3 .
Figure 3. Critical Rayleigh number as a function of fault width (without tapering) for three different fault heights, based on the theoretical analysis by Malkovsky and Magri (2016).

Figure 4 .
Figure 4. Lateral fluid pressure variation in the benchmark experiment (without tectonic regime consideration) is shown at a depth of 2.0 km The fluid pressure variation is compared in Figures 5-7 to the fluid pressure variation after the implementation of compressional, extensional, and strike-slip tectonic regimes, respectively.

Figure 5 .
Figure 5. 3D model results in a compressional tectonic regime (left) and associated fluid pressure variation (right).On the left, red colors correspond to positive temperature anomalies, and blue colors represent negative temperature anomalies.Fluid velocities are indicated, with red representing the fastest velocity and blue representing the slowest.On the right the fluid pressure variation, symbolized by a solid line, corresponds to the tectonic regime tested (here, compression).The fluid pressure variation symbolized by crosses corresponds to the pressure variation of the benchmark experiment.

Figure 6 .
Figure 6.3D model results in an extensional tectonic regime (left) and associated fluid pressure variation (right).On the left, red colors correspond to positive temperature anomalies, and blue colors represent negative temperature anomalies.On the right fluid velocities are indicated, with red representing the fastest velocity and blue representing the slowest.The fluid pressure variation, symbolized by a solid line, corresponds to the tectonic regime tested (here, extension).The fluid pressure variation, symbolized by crosses, corresponds to the pressure variation of the benchmark experiment.

Figure 7 .
Figure7.3D model results in a strike-slip tectonic regime (left) and associated fluid pressure variation (right).On the left, red colors correspond to positive temperature anomalies, while blue colors represent negative temperature anomalies.Fluid velocities are indicated, with red representing the fastest velocity and blue representing the slowest.On the right, the fluid pressure variation, symbolized by a solid line, corresponds to the tectonic regime tested (here, strike-slip).The fluid pressure variation, symbolized by crosses, corresponds to the pressure variation of the benchmark experiment.

Figure A1 .
Figure A1.Variation in permeability (x,y) incorporated within the fault zone of varying thickness.This heterogeneous variation in permeability (Equation3) corresponds to lithological variations (i.e., fault core/damage zone/basement) in the "multi-core fault" conceptual model(Faulkner et al., 2003(Faulkner et al., , 2010)).The chosen permeability values align with laboratory measurements conducted according to the Heap and Kennedy protocol (2016) on samples collected from the Pontgibaud crustal fault zone in the French Massif Central(Duwiquet et al., 2021b).Furthermore, these values are in line with the permeability data found in the Scibek (2020) database.

Table 3
Horizontal-To-Stress Ratio Values for Each Model