Abstract
The Pontgibaud crustal fault zone (CFZ) in the French Massif Central provides an opportunity to evaluate the high-temperature geothermal potential of these naturally permeable zones. Previous 2D modeling of heat and mass transfer in a fault zone highlighted that a subvertical CFZ concentrates the highest temperature anomalies at shallow depths. By comparing the results of these large-scale 2D numerical models with field data, the depth of the 150°C isotherm was estimated to be at a depth of 2.5 km. However, these results did not consider 3D effects and interactions between fluids, deformation, and temperature. Here, field measurements are used to control the 3D geometry of the geological structures. New 2D (thin-section) and 3D (X-ray microtomography) observations point to a well-defined spatial propagation of fractures and voids, exhibiting the same fracture architecture at different scales (2.5 μm to 2 mm). Moreover, new measurements on porosity and permeability confirm that the highly fractured and altered samples are characterized by large permeability values, one of them reaching 10-12 m2. Based on a thermoporoelastic hypothesis, a preliminary 3D THM numerical model is presented. A first parametric study highlights the role of permeability, stress direction, and intensity on fluid flow. In particular, three different convective patterns have been identified (finger-like, blob-like, and double-like convective patterns). The results suggest that vertical deformation zones oriented at 30 and 70° with respect to the maximum horizontal stress direction would correspond to the potential target for high-temperature anomalies. Finally, a large-scale 3D numerical model of the Pontgibaud CFZ, based on THM coupling and the comparison with field data (temperature, heat flux, and electrical resistivity), allows us to explore the spatial geometry of the 150°C isotherm. Although simplified hypotheses have been used, 3D field data have been reproduced.
1. Introduction
As potential geothermal reservoirs, crustal fault zones (CFZ) are still largely unexplored and unexploited (e.g., South Erzgebirge deep fault zone, Germany [1]; Gravberg, Sweden [2, 3]; Badenweiler-Lenzkirch Sutur, Germany [4]; Dewey’s Fault Zone, Canada [5]; Longmen Shan Fault Zone, China [6, 7]; Pasmajärvi Fault Zone, Finland [8]; and Luchlompal fault, Russia [9]). The cost-effectiveness of a reservoir depends on permeability, heat, and fluid flow. These conditions are often related to a specific geodynamic context and the physical conditions of heat transport.
A CFZ is generally not a simple plane but rather has a complex geometry defined by a succession of fault and fracture intersections reflecting intense deformation [10–12]. CFZ are manifestations of strain localization and modify the mechanical properties of the crust above the Brittle-Ductile Transition (BDT). These heterogeneities, characterized by the presence of a succession of undeformed protolith, damage zones, and fault cores, will cause a significant variation in the permeability of the crust [13]. Permeability variation in the crustal domain has been widely discussed since the 1990s [14, 15]. Overall, the permeability of a CFZ is controlled by (i) 3D distribution of structural features such as brittle planes or lithological interfaces, (ii) the physical properties (e.g., Young’s modulus and porosity) of lithologies along and adjacent to the fault zone, (iii) the state of regional stress, and (iv) more generally by the coupling between thermal (T), hydraulic (H), mechanical (M), and chemical (C) processes [10, 16–24]. At depths below the BDT, the deformation of granite leads to reductions in matrix porosity and, presumably, matrix permeability (e.g., [25]); however, these authors also note that microcracks can remain open in granite at high temperature (900°C) and high pressure (100 MPa), resulting in matrix permeabilities higher than that of the intact material.
Convective heat transfer is often invoked to explain temperature anomalies within faults [26–33]. Many parameters seem to condition the convective regime: the permeability ratio between the fault and its host, the CFZ permeability distribution, the fluid properties (viscosity and density variations), the CFZ thickness, and the depth and the dip of the CFZ [32–36]. The stress field is well known to influence CFZ permeability [20, 37–39] especially by influencing the void properties (fractures, pores, channels, and cracks). In fact, these voids under stress are expected to evolve over time and, according to several authors, favour shapes that will facilitate the circulation of fluids and stress propagation at different scales [40–45].
In order to estimate the geothermal potential of CFZ, it is necessary to consider a multidisciplinary approach [13, 15, 46–50]. Using such an approach (structural field measurements, laboratory porosity and permeability analysis, and numerical modeling) Duwiquet et al. [51] studied the geothermal potential of the Pontgibaud CFZ in the French Massif Central (FMC) and highlighted (i) the heterogeneous 2D distribution of permeability within the fault zone, (ii) CFZ with a near vertical dip as potential targets for high-temperature geothermal energy, and (iii) the position of the 150°C isotherm at a depth of 2.5 km. Although all these results were in agreement with previously acquired field data (geology, geophysics, and temperature measurements) [52, 53], the understanding of the 2D distribution of permeability has been made on a limited number of samples, and the modeling estimates have been made in 2D and based on TH coupling only. In addition, the structural framework poorly considered neighbouring faults to the Pontgibaud fault system. Using more field observations and more laboratory data, we aim to better characterize the geothermal system potential of the Pontgibaud CFZ in order to establish prospecting guides for geothermal exploration of CFZ worldwide.
Finally, to account for the mechanical aspect, we have performed a preliminary THM model where hypotheses are simple in order to focus on one single but important issue: could mechanical stresses affect the hydrothermal convective regimes? The objective of this simplified 3D THM model is twofold: first, a simple synthetic model is dedicated to the demonstration of the effect of stress intensity and direction on fluid flow and thermal regime within the fault zone. Second, we implement mechanical effects within a 3D TH model of the Pontgibaud CFZ, with the unique objective of questioning consistency of the obtained results with field data. To our knowledge, THM models are generally applied to purely pressure-driven flow, where fluid density is constant and thus with no buoyancy-driven components (see, however, [54, 55]). However, Vallier et al. [24] used a similar approach as ours but with another objective being the study of mechanical effects on gravity anomalies, so that no information on the coupling between fluid flow and mechanical effects can be retrieved.
Thus, in this paper, we (i) provide additional qualitative and quantitative data necessary for the understanding of a hydrothermal system such as the Pontgibaud CFZ; (ii) insert these new data in a 3D geological model; (iii) explore by numerical modeling the effect of stress intensity, stress direction on fluid flow, and temperature anomalies; and (iv) evaluate in 3D the geothermal potential of the Pontgibaud CFZ.
1.1. Geological Setting
1.1.1. Geological Formation
Bounded to the east and west by the Aigueperse Saint Sauves fault and the Sillon Houiller fault, the “La Sioule” licence (Figure 1) is largely affected by late Variscan tectonics. Metamorphic formations and the Gelles granite, which are intrusive, characterize the field. Magmatic and hydrothermal veins are widely represented.
Covering a large part of the mapped region, metamorphic formations (green in Figure 1) can be separated into different lithologies. Biotite and sillimanite paragneiss extends southwards near the Miouze river and northeastwards to the Villelongue region. Orthogneiss associated with biotite and sillimanite paragneiss are present sporadically. Cordierite gneisses are more widely represented, but it is difficult to assign a cartographic boundary. Rare amphibolitic levels are interspersed in the cordierite gneisses. As an intrusive in metamorphic series, the porphyroid granite of Gelles (red in Figure 1) extends southwest of Pontgibaud, between the Sioule Valley to the East and the village of Tortebesse to the West. To the North, the last granitic points appear immediately to the west of La Goutelle.
The presence of thermomineral springs (such as Pranal and Ceyssat) is an indicator of present-day fluid circulation (Figure 1). Temperature measurements have been carried out in the Pontgibaud borehole between a depth of 55 and 120 m (International Heat Flow Commission database, http://ihfc-iugg.org/; borehole number FR198). From these measurements, a geothermal gradient of 41°C/km-1 has been inferred, and a heat flow of 110 mW·m-2 was estimated. Magnetotelluric campaigns have been conducted in the Pontgibaud fault zone area by TLS-Geothermics and IMAGIR between 2015 and 2017 [50, 51]. Two anomalies located in the brittle crust at depths of 3 and 8 km could be linked to the presence of magma, clays (smectites), and/or hydrothermal fluids. These two anomalies appear to be related to the Pontgibaud fault zone.
1.1.2. Structural Setting
In order to understand the spatial organization of potential drains of some geological structure, it is necessary to first consider the tectonic evolution of the study area (Figure 2). We will then investigate the effects of the direction of the principal stress on the fluids flow in preexisting fractures. (1)Famennian to late Namurian (365 Ma-315 Ma) From Famennian to the early Visean, the maximum main horizontal stress axis is rotated 40° clockwise from NNW-SSE to NNW-SSW. During the Visean, the maximum main stress axis spatially varied from the vertical (extension, detachment of Argentat) to the horizontal (strike-slip, sinistral kinematic of the “Sillon Houiller”). Then, during the late Namurian, the maximum main horizontal stress shortening axis continued to rotate clockwise until it reached a NE-SW direction (Figure 2).(2)Early to late Stephanian (305-295 Ma) During the lower and middle Stephanian, N-S compression created a system of conjugated N140°-170°E sinistral and N10°-50°E dextral with N80°120°E reverse faults. This deformation episode was responsible for the NE-SW sinistral fault of the “Sillon Houiller.” During the transition from the middle to the upper Stephanian, the axis of rotation of the compression changed from N-S to NW-SE (Figure 2). This rotation activated or reactivated dextral conjugated faults N 120°E and sinistral N160-180°E. During the late Stephanian, the compression became E-W and resulted in a strong folding of coal deposits along the Sillon Houiller fault [59].(3)Permian to Jurassic (295-205 Ma) Following Stephanian compression, the N-S/NE-SW Permo-Triassic extension is represented in the basement by normal faults, tension cracks, and veins [59–62]. On the Cévennes border, this N-S extension followed NW-SE directions inherited from the Stephanian, forming NE faults such as the Cévennes fault [59, 63]. At the Trias-Lias border, a major metallogenic event occurred in the Massif Central, with the formation of a large number of veins and stratiform deposits, whose age was determined at about 194 Ma [64].(4)Jurassic to Cretaceous (205-96 Ma) During the transition from the Triassic to the Jurassic, a NE-SW extension occurred in the Rodez and in the western Cévennes, creating normal synsedimentary faults [59]. The beginning of the NW-SE extension was accompanied by significant metallogenic and thermal phenomena. Studies show that the appearance of N-S- to E-W-oriented barite and fluorite veins in Morvan, Quercy, and Lozère can be attributed to the Early Jurassic [61, 66, 67]. The Jurassic-Cretaceous period corresponds mainly to the development of the Liguro-Provençal Ocean and then the Atlantic Ocean. The beginning of the Jurassic is characterized by a radical change in the directions of extension, changing from a period marked by a strong late extension to an extensive period where the stresses are much lower, causing small deformations. This N-S oriented extension episode was responsible for the appearance of normal N80°-130°-oriented faults. The rapid rotation of the direction of the minimum σ3 and intermediate σ2 stresses around the maximum subvertical stress σ1 led to a new N-S-oriented paleogeography that persisted in the early and mid-Jurassic (Figure 2). At the beginning of the Late Jurassic, the Dogger structure was destroyed by another rotation of the minimum σ3 and intermediate σ2 stresses around the maximum vertical stress σ1. The previous E-W extension disrupted the organized N-S paleogeography elevations of the Middle Jurassic.(5)Cretaceous to late Pliocene (96-2.5 Ma) The establishment of a rift characterizes the Cenozoic in Europe: the European Cenozoic Rift System (ECRIS). This large rift system (about 1000 km) is associated with a series of transect faults throughout the Variscan massif. This produced a radical change in the stress regime throughout the Massif Central (Figure 2). After a long period of Permian-Mesozoic extension, a N-S compressive tectonic regime was established, the effects of which were visible from England to Bohemia [68–71]. In the Massif Central and its borders, the Late Eocene-Oligocene extension generated basins, grabens, and half-grabens bounded by normal faults, in which several hundred to thousand metres of continental sediments were deposited. The faults limiting these basins are generally N-S oriented but some NE-SW and NW-SE trends also exist. While the extension affected the entire Massif Central and its border, the deformation was concentrated near the regional fractures inherited from the Variscan orogeny. During the Oligocene, the Western European plate was affected by a system of N-S-trending grabens forming the Western European Rift. This continental rift can be divided into several sections linked by transforming regional faults with a NE direction [71, 72]. Among the N-S grabens, the Limagne half-graben, 15 km east of the PFZ, hosts a temperature anomaly that reaches 160°C at a depth of 3 km [73].(6)Present-day stress regime Several seismological stress-state evaluations were performed around the studied area. Mazabraud et al. [74] have determined five focal mechanisms from regional records below the Pulvérières and St. Gervais granites, located 20 km north of the PFZ. They exhibit reverse, normal, and strike-slip mechanisms. The maximum horizontal stress axis is oriented close to NW-SE. More recently, local stress states were evaluated during the exploration of the geothermal licence based on the seismological records acquired between 2015 and 2016 by TLS-Geothermics using 11 triaxial broadband seismometers CMG40 (see [53] for more details about this geophysical acquisition). Data were processed in 2018 by MAGNITUDE (a Baker Hughes company) to (i) identify microearthquake events, (ii) perform a focal mechanism analysis, and (iii) evaluate the stress state. 25 microearthquake events were identified and localized (mainly to the north of the studied area), and 14 focal mechanisms were determined (between -2 and -12 km asl). Confidence is very good for six Focal Mechanism Solutions (FMS) (good amplitudes fit at each station and focal spheres are well constrained), is good for three FMS, and is more uncertain for five FMS. They exhibit reverse (three FMS), normal (two FMS), and strike-slip regimes (nine FMS). The strike-slip regime occurred for all FMS with a high degree of confidence, for 2/3 of FMS with a good degree of confidence, and for 1/5 of FMS with a low degree of confidence. The nature of the stress tensor for each FMS quality set (see [74]) was then estimate by assuming that the following hypotheses are true: (i) the stress field is uniform within the studied area and (ii) the tangential traction on the fault plane is parallel to the slip direction. In each of the three cases, they obtain a strike-slip regime (i.e., with a subvertical intermediary stress axis and subhorizontal maximum and minimum stress axes). Directions of the maximum stress axis are bracketed between N135° and N150° for the high quality set and between N120° and N135° for the other two cases (Figure 2).
2. Materials and Methods
The multidisciplinary approach used here corresponds to a method widely used in oil and uranium exploration, CO2 capture and storage studies, and reservoir engineering. It combines a field approach (structural measurements and sampling of the Peyrouses 1 borehole), a laboratory approach (2D and 3D observations and measurements of permeability and connected porosity), and, finally, numerical modelling. The goal of our study is to integrate these field and laboratory measurements and observations into a 3D numerical model that is based on a Thermal, Hydraulic, and Mechanical (THM) coupling method using Comsol Multiphysics™.
2.1. Structural Observations
The structural measurements, performed in the field, make it possible to better characterize the intensity of the deformation in and around the Pontgibaud fault zone and also to put the geometry and architecture flow pathways into a 3D perspective. Sampling of the Peyrouses 1 borehole located within the Pontgibaud fault system provides a better understanding of the spatial distribution of pores and permeable fractures in this hydrothermal system, through thin-section observations and X-ray microtomography.
2.2. Laboratory Observations and Measurements
The samples used for laboratory measurements and observations were recovered from the Peyrouses 1 borehole (Figure 3), located within the PFZ (see Figure 1). These samples represent different facies: altered, fractured, altered and fractured, and intact facies (anatexite and microgranite). Thin-section observations on samples with different degrees of fracturing and alteration and at different scales help to characterize the petrology of the Pontgibaud hydrothermal system. 3D observations using a Phoenix nanotom microtomograph have been used at the Institut des Sciences de la Terre d’Orléans (ISTO). X-ray microtomography is a high-resolution, nondestructive 3D imaging technique. The samples were placed between a tungsten filament and a detector. During acquisition, the sample rotated 360° and the tungsten filament emitted X-rays onto the sample (using a voltage and current of 110 kV and 85 μA, respectively). The detector measured the attenuation of the emitted signal. A series of 2D images were produced, with pixel values dependent on the intensity of the emitted radiation and on the attenuation recorded by the detector. This attenuation is a function of the density and atomic number of the material through which the rays pass [76]. A 3D reconstruction of this object has been done, and each voxel reflects the average attenuation of the traversed material. However, the density contrast between two mineral phases may bias the information. A Beam-Hardening Correction (BHC) limits the loss of this information. After the reconstruction, the 3D volume was processed with the Volume Graphics VG Studio Max visualization software. This imaging technique provides high-resolution images (2.5 to 10 μm) and allows us to understand how fluids can flow in these heterogeneous materials. X-ray microtomography is a common method used to characterize oil [77], uranium [78], and geothermal reservoirs [79].
In order to quantify the heterogeneity of the system, permeability measurements have been performed at the Institut de Physique du Globe de Strasbourg (IPGS) in Strasbourg (France). For this purpose, 11 cylindrical samples (with a diameter of 25 mm and length between 25 and 40 mm) were prepared from an intact sample, an altered sample, a fractured sample, and an altered and fractured sample. The crystal size of these samples was small when compared to their diameter, as is recommended for samples used for measurements of permeability [79]. These samples were first oven dried under a vacuum at 40°C. Their connected porosities were calculated using the connected volume measured by a helium pycnometer (Micromeritics AccuPyc II 1340) and the bulk sample volume measured using digital calipers. Gas permeability was measured using a benchtop nitrogen parameter (see the schematic diagram in [81, 82]). All measurements were performed under a confining pressure of 1 MPa (samples were first left for 1 h at this pressure for microstructural equilibration) using the steady-state method. Volumetric flow rates were measured (using a Bronkhorst gas flowmeter) for six pressure differentials (measured using a Keller pressure transducer) to calculate permeability using Darcy’s law and to check whether the Forchheimer and/or Klinkenberg corrections are required (for more details on this protocol, see [83, 84]). The 3D numerical modeling discussed below will use these laboratory data.
2.3. Numerical Approach
Considering a very simple case, we will first try to understand the effects of stress direction and stress intensity and the effects of permeability on fluid flow, intensity, and depth of temperature anomalies. Then, we perform a simple preliminary study on the natural case of the Pontgibaud CFZ in order to get consistent results with field data. We first describe below the governing equations together with physical parameters. Then, the idealized simple model is described before the 3D geological model.
The poroelasticity interface of Comsol Multiphysics™ combines the regular fluid flow formulation with linear elastic solid mechanics of the porous medium matrix. Poroelasticity coupling means that the deformation of the porous medium affects fluid flow, while changes in volumetric stress will in turn affect the amount of motion and heat transport. This interface includes an expression of the stress tensor as a function of the strain tensor and the Biot-Willis coefficient [85].
Mixing mass conservation of the fluid and momentum balance yields that the poroelastic model is governed by the following relation: where is the fluid density (kg·m-3), is the poroelastic storage coefficients (Pa-1, see Equation (7)), is the fluid pressure (Pa), t is time, is the fluid source term (kg·m-3·s-1), is the Biot-Willis coefficient (-), and is the volumetric strain of the porous matrix (-); Darcy’s law sets that is equal to
where is the permeability of the medium (m2), is the fluid dynamic viscosity (Pa·s), and is acceleration due to gravity (m·s-2).
The governing equation for the momentum balance model is
Here, (Pa) is the stress tensor and (kg·m-3) is the total density. Notice that negative compression convention is considered.
The components of the strain tensor depend on the displacement vector (m). Under the small strain assumption, we have
The stress-strain relationship for linear materials relates the stress tensor and the strain tensor through Hooke’s law. Accounting for the hydromechanical coupling, we have for isotropic materials. Here, is the second order identity tensor (-), is the volumetric strain (-), and and are Lamé’s parameters (Pa) and are a function of Young’s modulus (Pa) and Poisson’s ratio (-) by the following relationship:
The storage coefficient (in Equation (1)) depends on the porosity (-), Biot coefficient (-), fluid compressibility fixed at [86], and matrix compressibility (Pa-1):
The matrix compressibility is fixed with Young’s modulus and Poisson’s coefficient by the following relationship:
Consequences of Equations (7) and (8) on numerical values of are the following: the coefficient is time-dependent through porosity variations and the obtained values vary between 10-11 and 10-10 Pa-1. These values are consistent with the compressibility values provided by Freeze and Cherry [87].
The poroelastic material model obeys Equation (3) which describes changes in the stress tensor and porous matrix deformation due to boundary conditions and changes in pore pressure. Equation (3) defines a state of static equilibrium because the changes in the solid equilibrate quickly, implying that there are no time-dependent terms. Still the time rate of change in strain appears as a coupling term in Equation (1) because the solid equation becomes transient when solved simultaneously with a time-dependent flow model.
Details of the coupling equations between Darcy’s law and Fourier’s law can be found in Guillou-Frottier et al. [88]. For our 3D experiments, realistic fluid and rock properties have been used. The fluid density law corresponds to fitting experimental data for pure water and for temperature and pressure ranges of 0-1000°C and 0-500 MPa, respectively (NIST Chemistry WebBook [89, 90]): where is the temperature (°C) and is the pressure (Pa). The dynamic viscosity of the fluid is a function of temperature , such that [91, 92]
Some rock and fluid properties are set as constant for the parametric study and in the large-scale numerical modeling of the Pontgibaud system. Their values are summarized in Tables 1, 2, and 3 (see Supplementary Material).
Permeability is a critical parameter in the functioning of geothermal systems and will be our main adjustment variable. We will distinguish between the permeability of the basement and the permeability of the fault zone, whose variations will be differentiated. In the basement, and up to the BDT, the permeability will be subject to the lithostatic pressure and will therefore be considered as decreasing with depth. It will be characterized by the following equation [93]:
Here, (m2) is the depth-dependent permeability of the basement and (m2) is the maximum permeability (present at the surface) taken here at . (m) is the depth below sea level, and (m) defines the decrease in permeability with depth, taken here at 2500 m, following the models of [51, 88, 93].
When considering the fault zone permeability , we have to consider the depth variation (as Equation (9)) but we must also consider the lateral variations observed and measured in the field and in the laboratory. This function of permeability is similar to that already used in [51], where the permeability varies in 2D along the and axes with a decrease in permeability with depth. In order to consider the permeability in 3D, we add a “” coordinate to this function and we then obtain where () is the permeability of the fault dependent on the three space directions and will be the maximum permeability at the surface. The choice of 100 and 101 means that the permeability will vary over two orders of magnitude. The permeability here will vary according to a sinusoid applied to the three directions . The term λ corresponds to the wavelength of the sinusoid. In order to approximate the observed and measured permeability variation (see below), we take a small value of λ.
The numerical simulations have been performed with Comsol Multiphysics™ in two different stages, illustrated in Figures 4 and 5, respectively. The first one is aimed at studying the effects of mechanics on fluid dynamics, through simplified and idealized hypotheses, while in the second one, previous geological observations and measurements are incorporated in a large-scale 3D THM model.
(a)
(b)
The first stage consists of investigating the effects of the fault direction, of the fault permeability, and of the stress value on the fluid flow, on the intensity of the temperature anomaly and on its depth. To do this, a series of numerical experiments have been performed with three variables: directions of the vertical deformation zone that vary from 0° to 90° (with a 10° step), stress values that vary from 0 to 100 MPa (with a 10 MPa step), and permeability values that vary from to (with a step) (Figure 4). The mesh consists of 175,789 tetrahedra, with a maximum size of 200 m for low permeability values and a maximum size of 20 m for high permeability values. The stress is applied on the red side (see Figure 4); the other three lateral borders have roller conditions (no perpendicular displacement). At the base of the model, an embedding condition (no displacement) and a heat flow of 80 mW·m-2 are applied. The vertical sides are thermally and hydraulically insulated, and at the surface, a temperature and pressure of 10°C and 105 Pa are imposed, respectively. A positive pressure evolution with depth is imposed as an initial condition.
The initial thermal conditions of these models are defined from the stationary state obtained for a pure conductive regime. At , values of permeability are assigned to the different compartments and the fluid flow is initiated.
The second stage concerns the numerical modeling of the Pontgibaud hydrothermal system in 3D (Figure 5). Here, the 3D geometry is constrained by the field observations and measurements. Indeed, based on the tectonic history (see Figure 2) and the measurements acquired in the field, different fault families have been individualized. The complete mesh consists of 244,475 tetrahedra, with a maximum size of 200 m and a maximum size of 10 m in most permeable zones. The mesh is refined for the highest permeable zones as well as for fault and lithology intersections. The result of this numerical modelling will be compared with field data, such as heat flux and geothermal gradient, as well as resistivity data from [53]. Considering the thermal relaxation time of magmatic chambers [94], calculations were performed up to a time of .
One of the critical points in setting up regionally scaled mechanical models are the boundary conditions [96, 97] that result from the far-field stresses. The overall stress orientation is often hinted from the integration of several measurements in a given region (such as well tests and GPS measurement). In our case, the present-day stress orientation in the surroundings of the Pontgibaud CFZ is approximately N145° (see section “structural setting,” above). As far as magnitudes are concerned, few or no indications usually exist, and our case is no exception. In most cases, numerical models include the boundary magnitudes as unknowns which are estimated through sensitivity analyses [98] or inversion studies [96]. In our case, we propose to take advantage of in situ measurements from the Chassoles borehole located at 60 km south-east of Pontgibaud ([95] Figure 4). This well was drilled for the “Géologie Profonde de la France” research program. The stress field has been determined by hydraulic testing. The well encountered only orthogneiss and intersected several subvertical fault zones at depths of 180 and 630 m. A total of 22 tests were run, nine between 340 and 470 m and 13 between 734 and 853 m. Many tests yielded fractures. At Chassoles, down to 7 km, the maximum stress is vertical but nearly equal to . This stress profile has been implemented as boundary conditions for the 3D model of the Pontgibaud CFZ. The heterogeneities in rock properties and structural features affect the stress distributions [99], and the extrapolation of measurements to other locations and/or depths appears to be a strong assumption. Still, regarding the high number of unknowns and parameters in our model, we decided to keep these observations as inputs in our first computations for the sake of simplicity.
3. Results
3.1. Structural Characterization
The north-south structures (Figure 6, A, C, D, and H) represent fairly clear topographic boundaries and underline the axis of Pb-Zn mineralized veins. The dips of these structures vary between 60 and 90°. The alignment of this axis is often interrupted by transverse structures and masked by the abundance of lavas in the eastern part. Although masked by the alluvium of the Sioule Valley, these north-south structures (Figure 6, C and D) underline a homogeneous alignment over 40 km.
E-W structures also mark the topography (Figure 6, E and K). With a dip between 60 and 90°, field observations and the geological map show an occurrence of these structures every 10 km in the southern part of the FMC and near the Sillon Houiller.
The NW-SE (Figure 6, G, I, and J) structures are quite frequent and often overlap the E-W faults without inducing significant offset.
The NE-SW (Figure 6, B, D, and F) structures can be seen widely on the geological map, from the west of the Sillon Houiller to Limagne. They are particularly well expressed on the Gelles granite. Mapping observations show that these faults are not too affected by N-S faults. The dip of these fault planes are mostly vertical.
To sum up, the Pontgibaud CFZ is a large deformation zone (>3 km) extending over 40 km, mainly characterized by north-south structures, and crossed by E-W, NE-SW, NW-SE, and some other NNE-SSW structures. These field observations have been used to determine the 3D geometry of our Pontgibaud large-scale numerical model. Before that, we will visualize and quantify the permeability distribution of the PFZ.
3.2. Thin-Section Observations (2D)
The following section summarizes the main lithologies present as well as the void space arrangement, for different scales of 2D observations. In these images, the voids are represented in white and the fractures in red. In order to make sure that the voids present are not caused by thin-section manufacture, we have highlighted the voids where secondary mineral phases have precipitated.
The facies of the samples are essentially fractured and altered (Figures 7(b), 7(c), 7(d), and 7(f)). Quartz, feldspars, biotites, and muscovites are mainly found. There are also secondary mineral phases (such as tourmaline, pyrite, arsenopyrite, scorodite, and sericite) and accessory minerals (such as apatite and zircon). The least altered and unfractured sample is shown in Figure 7(a) (LPA view). Quartz are most often xenomorphic and can form clusters of small crystals or develop over a wide range with an automorphic tendency, giving the rock a porphyritic texture. Biotite is also observed and generally automorphic and shows an onset of chloritization. We also observe the polysynthetic macle of plagioclase in Figure 7(a) (LPA view).
(a)
(b)
(c)
(d)
(e)
(f)
Within these altered and fractured facies, voids are present within cores of quartz mineralization (Figure 7(b)) but also at the walls of these mineralizations (Figures 7(b) and 7(c)). While some fractures have voids through which fluids can flow (Figure 6, D and E), other fractures are completely filled with clays, secondary mineral phases, and/or oxides. In addition, fractures can also be cemented by cataclasites (Figure 7(f)). These cataclasites consist of angular quartz clasts, plagioclase, sericite clays, and some calcites. Due to their similar petrological properties, most of these angular clasts should originate from the basement.
Most voids and secondary mineral phases appear to be associated with the development of fractures and cracks. Their interconnection could serve as potential drains for fluid circulation. The architecture of this interconnection at different scales could have a key role on the ability of the system to allow fluids to circulate.
Figures 8(a)–8(d) show images of the same sample but at progressively smaller scales (from 2 mm to 50 μm). These zoomed images are focused at the core of a fracture network in order to understand the spatial distribution of voids in the central part of the fracture network. Figures 8(e)–8(h) again represent the same sample but at progressively smaller scales (from 500 μm to 50 μm). These zoomed images are focused at the extremities of a fracture network in order to see how fractures and pores spread to the extremities of these networks.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
At a scale of 2 mm, Figure 8(a) shows macroscopic voids. They are distributed within the quartz vein and also at the walls of this mineralization. The secondary minerals and other smaller voids appear to form a network of fractures (red line) that extend away from the quartz vein, forming a dendritic network. At a scale of 500 μm (Figure 8(b)), other voids can be seen, still marked by secondary mineralization, which also follows the dendritic fracture network. By increasing the zoom again, we observe voids (Figure 8(c)) that we could not distinguish in Figure 8(b). These voids are still present within this dendritic network and are still marked by secondary mineralization. By increasing the zoom again, up to 50 μm, we still notice voids that we could not distinguish at the previous scale. These voids also follow a dendritic pattern. We find an interconnection of fractures at different scales, which form a dendritic fracture network. The precipitation of secondary mineral phases inside this dendritic fracture network shows that fluids have circulated.
Figure 8(e) shows another sample with a fracture network organization identical to Figure 8(a). Fractures marked by secondary mineralization, or by the presence of voids, appear to propagate spatially. By zooming in on Figure 8(e), we find voids that we could not distinguish previously (Figure 8(f)). These voids seem to follow the fractures and are also marked by the presence of late mineralization. Zooming in again on Figure 8(h), we now see a distribution of fractures that was not visible in Figure 8(g) and with an architecture not yet observed. This propagation of fractures across the box is in the form of rectangular drainage patterns (Figure 8(h)).
These flow patterns are widely present in natural drainage systems. These observations give us an overview of the spatial representation of the fractures and voids. We will now see if these structures are connected in 3D.
3.3. X-Ray Microtomography Observations (3D)
The visualization of fractures and 3D voids in a sample taken from the Pontgibaud CFZ (see location in Figure 1) are shown in Figures 9(a)–9(d). The sample is characterized by a fractured part as well as a matrix part. Samples were cored at the edge between the fractures and the matrix (Figure 9(b)) as well as in the matrix (Figures 9(c) and 9(d)).
The fractures are shown in red (Figure 9(b)), and the voids are shown in yellow (Figures 9(b)–9(d)). The fractures extend over the entire length of the sample and are heterogeneously distributed (Figures 9(b)–9(d)) in space. With a resolution of 10 μm, this 3D illustration of a complex architecture of the fracture network shows an interconnection of the voids present in these fractures (Figure 9(b)) that could facilitate the circulation of fluids. The 3D visualization of the matrix (Figure 9(c)) highlights the voids with a resolution of 10 μm. At this resolution, one observes globular clusters that appear to be interconnected. With a finer resolution of 2.5 μm, we can see more pores present in the matrix (Figure 9(d)). This resolution allows us to confirm the presence of voids that form a planar structure.
Thin-section (2D) and X-ray microtomography (3D) observations allow us to observe the architecture of the fracture networks and the distribution of pores at different scales both within the fracture networks and in the matrix. The porous flow structures of the Pontgibaud hydrothermal system have multiple scales and a nonuniform distribution of the space available for fluid circulation. The anisotropic nature of the distribution of permeability in the different dimensions of space has thus been quantified by measurements of permeabilities and connected porosities.
3.4. Porosity and Permeability Measurements
The data show that there is a general trend of increasing permeability as a function of increasing porosity (Figure 10). Connected porosity values range from 6 to 22% (Figure 10). In most cases, the altered samples (with and without fractures) have the highest connected porosity values (around 20%), and the intact samples have the lowest porosities (between 6 and 8%). Permeability values range from to . Altered samples with fractures have the highest permeability values (between 10-13 and ). Samples that are either altered or contain a fracture are characterized by intermediate permeabilities, between 10-16 and . The intact samples have the lowest permeability anisotropy, which is most pronounced in the samples that contain fractures (with and without alteration) (Figure 10). For example, the unaltered sample containing a fracture (Figure 10) is a factor of 20 more permeable in the direction (blue dot in Figure 10) than in the direction (green dot in Figure 10) (despite their similar porosities). The permeability of the altered sample containing a fracture (Figure 10) is three orders of magnitude more permeable in the direction (green cross in Figure 10) than in the direction (blue cross in Figure 10), despite their similar porosities.
3.5. Numerical Approach
More than 2500 numerical experiments were run to understand the effects of fault direction, stress intensity, and permeability on the intensity and depth of the positive temperature anomaly. Figure 11 synthetizes these results. This figure presents the value of the positive temperature anomalies, as well as their depths, as a function of the stress intensity (horizontal axis), the angle with respect to the stress direction (vertical axis), and the permeability ratio between the fault and the host rocks defined as
(a)
(b)
(c)
This permeability ratio compares the maximum value of the permeability of the fault with that of the basement (here fixed at ). Other . values were tested, but higher values do not correspond to the field data, and lower values of permeability do not change the results.
3.5.1. Numerical Parametric Study
By considering all the directions of the deformation zone and all the stress values, we can see that the largest temperature anomalies and the smallest depths will be concentrated for values of [50-200] and [200-400] (Figure 11(a)). By varying the direction of the deformation zone with respect to the stress, for permeability ratios [50-400], the higher temperature anomalies and the smallest depths will be concentrated in directions of 30 and 70°. This tendency is, however, less visible for [1-50]. By varying the stress, we notice that for [1-50], the intensity of the stress does not significantly influence the value of the positive temperature anomaly as well as its depth. However, trends seem to emerge for stresses of 30 and 90 MPa. These stress values seem to concentrate the highest temperature anomalies for the directions of 30 and 70° relative to the stress as well as for [200-400]. For values of [50-200], we do not observe significant effects of stress on the depth and intensity of the temperature anomaly.
Finally, for [50-400], the deformation zone directions of 30 and 70° relative to stress will concentrate the highest temperature anomalies at the shallowest depths. For these directions, stresses of 30-90 MPa will have the highest temperature anomalies.
Figures 11(b) and 11(c) show the fluid flow for [50-200]. For Figure 11(b), the direction of the deformation zone is oriented at 30° to the stress, which is 30 MPa. The fluid flow (marked by the black arrows) is defined by an upward movement in the center and downward movements at both ends. The bulging of the 70°C isotherm (shown in red in Figure 11(b)) is the consequence of this fluid circulation, which follows a “blob-like” pattern. For Figure 11(c), the fault direction is 70° with respect to the stress, which is 50 MPa. Fluid flow is also defined by an upward movement in the center and downward movements at both ends. However, compared to Figure 11(b), the downward movements are more concentrated around the upward movement. The morphology of the 70°C isotherm (shown in red in Figure 11(c)) is the consequence of this fluid circulation which follows a “finger-like” pattern. As can be seen, there seems to be different types of convection patterns for different stress and fault direction values. Note that the irregular morphology of the red surfaces may be due either to the 3D heterogeneity of the permeability or to the effects of stress direction and/or intensity. To solve for this important issue, convective patterns have been systematically observed for all experiments.
Observations of the convective figures from the numerical simulations make it possible to distinguish three zones as a function of the stress intensity, fault direction, and the value of the temperature anomaly (Figure 12). The brown zone and the green zone correspond to a ratio [1-200]. The red zone corresponds to a ratio [200-400]. For this permeability ratio, the temperature anomalies will be between 88 and 150°C. For the lowest permeability ratio ( [1-200]), temperature anomalies are between 7 and 88°C. For these values, two convective regimes are observable. Firstly, for stress values below 40 MPa and for a fault direction below 50°, the blob-like pattern described in Figure 11 is found. This convective regime occurs along two dimensions of the space (along the and axes). For stresses greater than 40 MPa and for a fault direction greater than 50°, the finger-like pattern described in Figure 11 is preferred. This convective regime is expressed along two dimensions of the space (along the and axes). For the highest permeability ratios ( [200-400]), we find a continuous upward and downward motion in all three dimensions of space. Circulation of warmest fluids is along the and axes, while the circulation of the coldest fluid is along the and axes. This circulation of fluids occurs according to a “double-like” pattern.
By considering a simple geometry and idealized hypotheses, it can be seen that the stress direction, the stress intensity, and the permeability ratio between the fault and the basement are parameters that have an influence on fluid flow and consequently on the intensity and depth of the temperature anomalies. Considering now a more complex geometry, constrained by field data, we can try to see if similar trends as those found previously (finger-like, blob-like, and double-like convective patterns) can be retrieved in a 3D THM numerical model.
3.5.2. Large-Scale Numerical Modeling Validation
The results of the numerical modeling show that for a maximum permeability imposed on the Pontgibaud deformation zone of , the geothermal gradient would be 41.5°C/km (Figure 13(a)) and the heat flux would be 110 mW·m-2 (Figure 13(b)), at the location (Figure 1) where the measured and estimated field data are 41°C/km and 110 mW·m-2. The resistivity profiles show a low resistivity anomaly at a depth of 1 km (Figure 13(b)) and a larger low resistivity anomaly at a depth of 3-4 km. If the origin of this resistivity may be discussed, the comparison between the resistivity profiles and the temperature anomalies shows that the temperature anomaly within the Pontgibaud deformation zone is between 2 and 3 km deep. In the middle part of the Pontgibaud deformation zone and at the intersections of the different fault families (N-S, NE-SW, NW-SE, and E-W), the 150°C isotherm (in red on Figure 13(a)) would occur at a depth of 2.3 km.
Permeability variation (Figure 13(c)) imposed on large-scale numerical simulations corresponds to observations made in the field, in 2D thin-section images, and in 3D X-ray microtomographic images.
4. Discussion
4.1. Origin and Effect of Heterogeneities
Thin-section (2D) and X-ray microtomography (3D) observations have revealed spatially propagating fracture networks in which fluids can flow through pores at different scales, from 2.5 μm to 2 mm (Figures 7, 8, and 9). This fracture propagation has been studied on mineralized systems (St. Salvy, France, and Jebel Aouam, Morocco) and in which tectonic comminution processes have been identified [101–103]. Microscopic scale pores have also been observed in the walls of some mineralization. The most frequent process in the brittle part of the crust is hydrofracturing. This process depends on the temporal variation of the fluid pressure [104–106], causing a “break up” which may be at the origin of some voids. This process has been studied on sites with F (El Hammam-Morocco), Ba (Dreislar, Germany), and Ag and Au (Creede, Colorado, USA) mineralization [102]. Fluid circulation can also cause the dissolution of mineral phases in the basement. Dissolution may be accompanied by a loss of volume and thus cause the formation of pores within the matrix [90]. While the consequences of these processes can be observed in 2D and 3D, other processes can cause the opening or closing of the voids (such as volume reduction, volume expansion, and pore collapse), thus influencing the capacity of the system to allow fluids to circulate (Figure 14).
(a)
(b)
(c)
Measurements of permeabilities and connected porosities from samples from Pontgibaud show ranges of values between and and 6 and 22%, respectively (Figure 10). These values are within known ranges of permeability for rocks in a basement context [1]. Our data also highlight the heterogeneities of the permeability of a hydrothermal system, which could be caused by the processes of tectonic comminution, hydrofracturing, and mineral phase dissolution, processes that can form structures and pores at different scales. The permeability of a rock mass will be very much controlled by the fracture density and, importantly for hydrothermal systems, the proportion of open fractures [82, 84, 107, 108]. Our permeability data also show that alteration can increase permeability (Figure 10). Previous studies on volcanic hydrothermal systems have also shown that hydrothermal alteration can greatly influence matrix permeability [109–113]. For these natural hydrothermal systems, permeability should not be approximated to a single value but should rather be represented as pipes or channels embedded in the neighbouring system [114].
According to Bejan and Lorente [115], architectures at different scales and arranged in a hierarchical pattern would promote fluid flow. These natural drainage systems (Figures 8 and 9) may have lower flow resistance than a similarly porous volume. The constructal theory [40] allows us to describe the natural tendency of flow systems to generate multiscale shapes and structures favouring fluid flow. It has also been shown that the transition between heat transfer at the nanoscale and more conventional scales is governed by the smallest dimensions and by their spatial heterogeneities [116]. The origin of these heterogeneities could be the stress (with tectonic comminution) and the water (with fluid pressure). Finally, with knowledge of current parameters, the geometry of the Pontgibaud CFZ might correspond to an optimal one, which facilitates efficiency of convective flow.
4.2. Effect of Stress Direction and Stress Intensity on the Fluid Flow
Multiple parameters can influence fluid flow, and the 3D behaviour of convective cells can be complex [117]. López and Smith [31] showed that temperature anomalies can be directly influenced by the permeability ratio between the fault and the host. 3D modeling based on TH coupling has shown the influence of this permeability ratio on the morphology of the temperature anomalies [33]. Our parametric study (Figures 11 and 12) made it possible to approach three factors influencing fluid flow: stress intensity, stress direction, and the permeability ratio between the fault and the basement. For high permeability ratios between [200-400], the poroelastic effects due to stress state variation have no significant effect on fluid flow pattern. However, for lower permeability ratios [1-200], changes in the fluid flow pattern can be observed depending on the stress direction and stress intensity.
4.3. Positive Temperature Anomaly on the Pontgibaud CFZ
As shown in Figure 13, the results of Pontgibaud large-scale numerical modeling compared with heat flux, geothermal gradient, and resistivity data show that the 150°C isotherm could be located at a depth of 2.3 km within the Pontgibaud CFZ. This upwelling of the 150°C isotherm is the result of fluid circulation through a blob-like convective pattern. This will produce a temperature anomaly within the Pontgibaud hydrothermal system, but lateral thermal diffusion is observed within the fault zone itself, thus defining a larger geothermal reservoir.
Indeed, considering a simplified 3D geometry, this lateral temperature diffusion has been highlighted in recent numerical models of fluid flow around faulted zones [33].
The Pontgibaud CFZ is located at 45° to the main regional stress (Figure 14). In this direction, the convective regime diagram (Figure 12) predicts blob-like convective pattern fluid flows for ratios [1-200], as well as for stress intensities below 40 MPa. In view of these parameters, the result of the convective regime diagram seems to be consistent with the validation of the large-scale numerical model of the Pontgibaud CFZ in which blob-like convective pattern fluid flows would be present.
If the consistency between general concepts from the very simplified numerical modeling is effectively validated by natural system results, then it would be possible to use the convective pattern diagram in order to characterize convection regimes of a natural system according to a tectonic regime.
Figures 13(a) and 14(c) show a positive temperature anomaly that extends widely on both sides of the Pontgibaud CFZ. The origin of this anomaly would be a blob-like convective pattern. Both blob-like and finger-like convective patterns are found for low permeability ratios (see Figure 12). Finger-like convective patterns are characterized by isotherms covering a smaller space than blob-like convective patterns. Consequently, the volume of the temperature anomaly will be less in a finger-like convective pattern than in a blob-like one. As we have seen, the parameters influencing the formation of these both convective pattern are the intensity and the direction of the stress. Thus, in a geothermal exploration context, in geological settings where the permeability is quite low, it would be advantageous to target areas where the intensity of the stress is quite low and where the direction of the stress is less than 50° with respect to the deformation zone. In this case, there will be more widespread temperature anomalies.
For high permeability zones, the intensity and direction of the stress are no longer parameters influencing the convective regime under thermoporoelastic conditions. The double-convective pattern gives higher temperature anomalies than the two patterns previously described (). Although the temperature anomaly is more intense, the deformation of the isotherms is more restricted than for blob-like convective patterns. Thus, in areas of high permeability, the temperature anomaly will be the most intense, but the volume covered will be less than for blob-like convective patterns.
4.4. Limitation and Perspectives
The theory of thermoporoelasticity considers the two effects of temperature change and fluid pressure on the elastic behaviour of rocks. This theory could be further treated in mechanical terms, in which slip and dilation tendency could be calculated [118, 119]. In fact, some fault segments subjected to critical stresses would act as conduits for the circulation of fluids [120]. As a consequence, some preferential flow paths can appear under the action of mechanics.
The synthetic model illustrates the impact of mechanics with one boundary stress only. The effects are expected to be greatly emphasized under anisotropic boundary conditions because of the shearing conditions they introduce.
Permeability is widely recognized as a central parameter in the understanding of a geothermal system. It is also a factor that evolves in space and time. Stresses, fluid overpressure, thermal cracking, mineral dissolution, and precipitation are just a handful of examples that can affect the permeability of a geothermal system [121–125]. All of these phenomena help to explain how a system that is initially impermeable evolves into a geothermal reservoir potential over time. Another method to improve such 3D models would be to be able to consider permeability as a dynamic parameter that varies in response to different geological processes, including varying tectonic regime as illustrated by the Pontgibaud CFZ ().
The parametric studies were realized by considering the faults as vertical, faults which from the 2D parametric study on TH coupling showed that the highest temperature anomalies concentrate at the lowest depth. However, considering a THM coupling, it would be interesting to be able to consider the dip of the structures. It would then be possible to highlight stress directions, permeability ratios, stress intensity, and dips as favourable targets for high temperature geothermal energy. In general, geothermal exploration must consider the intensity of the stress, the direction of the stress, and the variability of permeability.
Finally, all these results are only applicable for large deformation zones in which temperature anomalies are induced by buoyancy-driven flow and not in zones where there is significant topography. Moreover, we only consider one fluid phase (liquid water). If we were to consider two phases, liquid and gas, the poroelastic effects on the compressibility of the gas could be an important process [126]. If the compressibility of the fluid can also influence the elastoplastic behaviour [127], the thermal dependence of compressibility might change the thermal distribution of the system [128].
5. Conclusions
Based on the results of this multidisciplinary study, we speculate that crustal fault zones represent potential reservoirs for high-temperature geothermal energy. Field measurements collected at the Pontgibaud crustal fault zone highlight a large N-S deformation zone, with a with a dip of 60° to the east and intersected by different families of NE-SW, NW-SE, and E-W faults. Laboratory observations show drainage structures with mainly dendritic geometry. The origin of this geometry here could be the fluid pressure and/or tectonic comminution. Such flow structures are widely present in different natural systems and have been suggested by the constructal theory as optimal geometries for fluid flow. All these observations and measurements have been inserted in 3D numerical models with THM coupling. The first simplified numerical study allowed us to conclude that (i)for the simplified used boundary conditions, vertical deformation zones, having an angle of 30° and 70° with respect to the boundary stress, will concentrate the highest temperature anomalies, at the lowest depths(ii)different convective patterns could exist: finger-like convective patterns, blob-like convective patterns, and double-like convective patterns. These three convective regimes depend on parameters such as the intensity and direction of the stress, the permeability ratio between the basement, and the deformation zone
In a second step, large-scale numerical simulations of the Pontgibaud crustal fault zone were performed. Recalling that no calibration was conducted, these computations were ran essentially as a demonstration of their ability to study natural phenomena while incorporating field observations. Interestingly, they already show consistent results with field data. These results suggest that the 150°C isotherm is at a depth of 2.3 km and could therefore represent an interesting target for high-temperature geothermal exploration. This tendency will be more thoroughly studied in future works by using more constrained models (with, e.g., sensitivity analyses). Also, these results need to be confirmed by drilling by the operator. If the drilling confirms these results, then it could be used on other study zones such as Badenweiler-Lenzkirch Sutur (Germany), Longmen Shan Fault Zone (China), and Pasmajärvi Fault Zone (Finland).
Data Availability
Sample position, structural measurements, permeability and connected porosity measurements are included in the paper (Figures 3, 6, and 10). The physical properties of fluids and rocks are detailed in Tables 1, 2, and 3 (see Supplementary Material).
Conflicts of Interest
The authors declare that there are no conflicts of interest for the publication of this article.
Acknowledgments
We sincerely thank the association “La Route des Mines Dômes-Combrailles” as well as the “Musée de la Mine” for allowing the sampling of the Peyrouses 1 borehole. We sincerely thank Laurent Bailly for his help on the thin-section observations and for his constructive remarks. We warmly thank Prof. Fabien Magri for these valuable suggestions and remarks. We are grateful to Séverine Caritg, Eric Fournier, Anne-Sophie Serrand, Antoine Armandine Les Landes, Simon Lopez, Arnold Blaisonneau, Bastien Hermant, and Jean-Michel Ars, for their fruitful discussions and relevant advice. The financial support was provided by TLS-Geothermics, ISTO, and BRGM and by the “GERESFAULT” project (ANR-19-CE05-0043).
Supplementary Materials
The supplementary materials contain three tables. They group together the physical and thermal properties of rocks and fluids inserted in the numerical modeling. (Supplementary Materials)