Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

On the geothermal potential of crustal fault zones: a case study from the Pontgibaud area (French Massif Central, France)

On the geothermal potential of crustal fault zones: a case study from the Pontgibaud area (French... hugo.duwiquet@gmail.com TLS‑ Geothermics, 92 The present study aims to understand the potential of a new and novel type of geo‑ Chemin de Gabardie, thermal play system for high temperature and electricity production: crustal fault zones 31200 Toulouse, France Full list of author information (CFZ). According to geological and geophysical data, the Pontgibaud fault zone (French is available at the end of the Massif Central) is suspected to host an active hydrothermal system at a depth of a few article kilometers. The deep geometry of the fault zone and the permeability distribution are the main unknown parameters that are required to assess the geothermal potential of the Pontgibaud site. Structural and thin‑section observations, laboratory permeability and connected porosity measurements and X‑ray micro ‑tomography observations suggest that the hydrothermal system behaves like a double matrix‑fracture perme ‑ ability reservoir. Numerical modeling in which we varied the fault dip and the ratio between the fault zone permeability and host rock, R, was performed. Results indicate that three main convective regimes can be identified (weak convection, single cellular ‑ type convection and bicellular convection). For a sufficiently high fault zone permeabil‑ −15 2 ity (> 1 × 10 m ), buoyancy‑ driven flow creates a positive thermal anomaly of several tens of °C at a depth of 2–5 km. For a vertical fault zone, the thermal anomaly is larger for higher R values. Numerical models, then applied to the geologically constrained Pontgibaud fault zone, show that a temperature of 150 °C at a depth of 2500 m can be −14 2 obtained for a fault zone permeability of 1.6 × 10 m . Based on a multi‑ disciplinary approach, this work establishes a potential predictive tool for future high‑temperature geothermal operations within basement rocks hosting large‑scale fault systems. Keywords: Crustal fault zones, Permeability, High‑temperature geothermal system, Numerical modeling, Pontgibaud, French‑Massif Central Introduction Exploring new geothermal targets requires the understanding of the factors affecting fluid circulation and heat transfer (Rowland and Sibson 2004; Fairley 2009). Understand - ing the distribution of permeability in the crust remains an essential component for the general comprehension of a geothermal model in the crustal domain, and therefore for the success of a geothermal prospect (Moeck 2014). Nowadays in France, most geother- mal exploration licences are located in sedimentary basin areas or within graben struc- tures. The Paris and Aquitaine basins, characterized by simple sedimentary structures, allow geothermal exploitation from low (30 °C) to medium (110 °C) temperature systems. © The Author(s) 2019. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 2 of 29 Geothermal exploitation is also rife in the northeast of France in the Upper Rhine Graben. The enhanced geothermal systems at Soultz-sous-Forêts and Rittershoffen (Alsace, France) reach temperatures up to about 160  °C (e.g., Pribnow and Schellschmidt 2000; Baujard et al. 2017). However, no exploitation is currently being developed in the basement domain that host crustal-scale faults that are likely to modify fluid circulation. Numerous studies highlight that crustal-scale fault zones located in an extensional regime represent efficient conduits for meteoric fluids to flow down to mid-crustal depths (Barton et al. 1995; Faulds et al. 2010; Mckenna and Blackwell 2004; Siebenaller et al. 2013; Haines et al. 2016; Roche et al. 2018a, b). At these depths (8–15 km), corresponding to, or close to, the brittle–duc- tile transition, temperatures may exceed 400 °C (Violay et al. 2017). Fault zones appear as areas of high permeability in which fluids can circulate and focus temperature anomalies at the surface (Faulds et al. 2010; Taillefer et al. 2018a, b). However, the generation of a geo- thermal resource in a fault zone likely depends on its geometry with respect to the present- day stress field, on its permeability structure and on the local thermal regime. The present study aims at understanding the potential of a new and novel type of geothermal system for electricity production: the crustal fault zone (CFZ). Crustal fault zones can be defined as crustal-scale heterogeneities, which localize the deformation and thus modify the mechanical properties of the crust, down to the brit- tle–ductile transition (BDT), whose depth depends on the local thermo-mechanical regime. In some cases, CFZ may cross BDT and extend below it (e.g., Famin et al. 2004; Jolivet et  al. 2004). CFZ should not be considered as single major fault but rather as a large deformation zone (several hundreds to thousands of meter wide). These crustal fault zones can be found worldwide and their geothermal potential can be highlighted by the occurrence of geothermal events, e.g., North Anatolian Fault Zone, Turkey (Süer et al. 2008), Liquine-Ofqui fault zone, Chile (Lahsen et al. 2010), Badenweiler-Lenzkirch Sutur, Germany (Brockamp et al. 2015). However, their geothermal potential has never been studied and even less quantified. Among the first studies on fluid circulation within a fault, Murphy (1979) proposed a theoretical approach using the Rayleigh number. He highlighted vertical thermal insta- bilities linked to convection cells. Convective heat transfer establishes positive and nega- tive thermal anomalies within a fault (Lowell 1975; Murphy 1979; Rabinowicz et al. 1998; Faulds et al. 2010). While many parameters seem to condition a convective regime, the permeability ratio between the fault and its host rock emerges as a key parameter (For- ster and Smith 1989; López and Smith 1995). While fault zones often host hydrothermal systems (Faulds et al. 2010), the permeability of these networks depends on a multitude of parameters. The relationships between (i) fault geometry and stress field (Barton et al. 1995); (ii) fractures and their connectivity, and (iii) alteration, dissolution and precipita- tion, are all processes that can cause faults to behave as conduits for (Caine et al. 1996; Mazurek 2000; Gudmundsson 2000; Grasby and Hutcheon 2001; Brogi and Fulignati, 2012; Belgrano et  al. 2016) or barriers to (Bense et  al. 2008; Bertrand 2016) fluid flow. Consequently, it is worth using a multidisciplinary approach to strongly constrain fluid flow models in hydrothermal systems (Oda 1986; Caine et  al. 1996; Bense and Person, 2006; Micarelli et  al. 2006; Faulkner et  al. 2010; Bense et  al. 2013). As an example, the recent study by Edel et al. (2018) predicts, using a combination of magnetic, gravity, seis- mic and geological data, a temperature of 150 °C at a depth of 2500 m in the northern Duwiquet et al. Geotherm Energy (2019) 7:33 Page 3 of 29 Vosges area (France). With an improved numerical modeling approach, Guillou-Frottier et al. (2013) reproduced and predicted temperatures at geothermal sites of Soultz-sous- Forêts and Rittershoffen. Since 2014, TLS-Geothermics, a French company involved in geothermal exploration has been keen to demonstrate the viability of CFZ as a geothermal exploration play for economic power generation. The demonstration of the potential of CFZ is underway at a site in the French Massif Central (FMC) using an exploration licence named “La Sioule”. Several geological observations and geophysical data, which are detailed below, represent good indicators for high-temperature geothermal exploitation (tempera- tures > 150 °C for shallower than 3500 m). The Pontgibaud fault zone, located in the La Sioule licence (Fig.  1) corresponds to a well-suited case study since numerous data are available (structural geology, lithology, topography, springs geochemistry, thermal properties, magnetotellurics, gravimetry, ambient noise tomography, and seismicity) (Bellanger 2017). Thanks to this complete dataset, and with new additional petrophysical data, 2D numerical modeling can help to better understand the links between fluid flow and permeability. The possible presence Fig. 1 Simplified geological map showing the “La Sioule” licence area (black contour) and a brief overview of the major stratigraphic units in the study area. French Massif Central (FMC). Coordinate system: WGS84 Pseudo‑Mercator EPSG:3857 (After the 1/50,000th geological map published by BRGM) Duwiquet et al. Geotherm Energy (2019) 7:33 Page 4 of 29 of convection cells is here examined to understand the establishment of positive or nega- tive thermal anomalies in the fault system. Permeability will be considered as a spatially variable parameter to analyze its effect on the hydrothermal regime of the fault system. Geological setting Geological formations The La Sioule licence area (black contour in Fig.  1) is located to the east of the Limagne graben in the Massif Central, in a Paleozoic basement domain and partially covered by the northern part of the Cenozoic volcanic formations (in blue in Fig.  1). It is bounded to the west by the Sillon Houiller fault zone, to the south by the volcanic massif of the Monts Dore, and to the north and east by the Carboniferous basin of Manzat and the granitic massif of Saint-Gervais, respectively. Variscan metamorphic formations (ortho and para-migmatites, orthogneiss, paragneiss and mica schists) are mainly present in the area (in green in Fig. 1). Carboniferous plutonic rocks are characterized by monzo- granite and granodiorite (red in Fig.  1) and Viseo-Serpukhovian volcano–sedimentary formations (brown in Fig.  1). Upper Pennsylvanian and Permian detritic sedimen- tary formation mark the late Variscan basins emplaced on the major crustal fault zone (mapped by the BRGM). Structures and favorable features for fluid circulation in the “La Sioule” licence The Sillon Houiller fault zone This fault zone, generally oriented N 20°, is located to the west of the area (Fig.  2) and is continuous at a scale of 400  km with a sinistral motion (Faure et  al. 1993, Joly et  al. 2008). Evidence in various volcanic deposits highlights the high permeability of this structure since the Permian. At both sides of the Sillon Houiller, P-wave velocity anoma- lies indicate that this structure affect the entire lithosphere (Granet et  al. 1995a). Such a large-scale structure, if permeable enough, would promote the ascent of hot hydro- thermal fluids from the base of crust to the surface. According to Sobolev et al. (1997), temperature varies between the different sides of the fault at 30  km depth, from 700 −2 to 925  °C. In addition, surface heat flow equals 113 mW  m near the Sillon Houiller (Fig. 2) (Lucazeau et al. 1984; Lucazeau and Vasseur 1989). This abnormally high value is thought to be related to a relatively high mantle heat flow (Lucazeau et al. 1984), which could be related to a supposed mantle plume (Granet et al. 1995b) or past subduction- related hot upwelling (Király et al. 2017; Roche et al. 2018ab). The Aigueperse–Saint‑Sauves fault zone This NE–SW fault system is more than 300 km long (Faure et al. 1993) and crosses the southeast corner of the licence (Fig.  2). It comprises a fault network, several hundred meters thick, and shifts the Sillon Houiller fault zone over a 31-km-long dextral-slip motion. The Puy de la Nugère and Puy Chopine monogenic volcanoes of the quater - nary Chaîne des Puys (CdP) volcanic field (located 10 km to the east of La Sioule licence area) are emplaced on the Aigueperse–Saint-Sauves fault zone. The presence of Ceyssat, Duwiquet et al. Geotherm Energy (2019) 7:33 Page 5 of 29 Fig. 2 Thermal, hydrothermal and geological features in the studied area. Location of Pranal and Ceyssat springs (brown dots), drilling of Peyrouses 1 (black dot), heat flow and geothermal gradient (red crosses, Lucazeau and Vasseur 1989; Vasseur et al. 1991; International Heat Flow Commission). Gelles granite (Négroni 1981). Heat production (Lucazeau 1981) and resistivity profile (black dashed line, Ars, 2018) Châtelguyon, Gimaux, Saint Myon or Montpensier springs characterizes this structure as a potentially open and efficient pathway for hydrothermal fluids. The Pontgibaud fault zone The N–S fault zone of Pontgibaud is parallel to the CdP. The CdP is a fissural volcanic field that hosts about eighty basaltic to trachytic volcanoes that erupted between 100 and 7 ka. This system (Fig.  2), more than 40 km long, is pinched between the Sillon Houiller and the Aigueperse–St-Sauve fault zones. In the studied area, it bounds the Viseo-Ser- pukhovian Gelles granite (Fernandez 1969; Négroni 1981). Some N–S faults of the Pon- tgibaud fault zone are filled by microgranite. Moreover, this N–S network is affected by As–Au, Pb–Ag and Ba–F mineralization, reflecting the ability of the system to transmit hydrothermal fluids during the Variscan. Some dating on mineralized bodies obtained Duwiquet et al. Geotherm Energy (2019) 7:33 Page 6 of 29 on phengite and illite suggest Permian to Jurassic ages (Bril et  al. 1991), which reflects late hydrothermal fluid circulation related to the opening of the Alpine Tethys and the Atlantic Ocean. This mineralization is attributed to a medium temperature phase (200– 400 °C) sub-contemporary to the end of the crystallization of the Gelles granite. If evi- dence for hydrothermal fluid circulation exists during the Carboniferous, but also from the Permian to Jurassic, current fluid circulation along this fault zone seems also very likely. Indeed, springs can be found along the present-day fault at several locations (see details below). Geothermal setting of the Pontgibaud area Springs and magma chambers Present-day fluid flow occurs within the Pontgibaud fault zone (Fig.  2). In addition, at the time of mining activities, the galleries were regularly flooded (Lebocey 2013). Due to cessation of mining activities in the 20th century, the mining wells were abandoned. Some wells, such as the one at Sainte-Barbe in Barbecot (4  km north of Pontgibaud), are marked by ion-rich and extremely CO -charged flows. In addition, the Pranal spring (Fig.  2) comes from a deep reservoir with estimated temperatures between 150 and 175 °C (estimated using Na/K geothermometer; Verney 2005). At a deeper level, thermo-barometric experiments on the CdP trachytes estimated the storage pressure and temperature conditions of magmas (Martel et  al. 2013). For temperatures ranging from 700 to 825  °C, the results of Martel et  al. (2013) indicate pressures between 3 and 3.5 kbar (i.e., between 11 and 13  km deep, respectively). The thermal relaxation time calculated for these reservoirs indicate that cooling started about 15,000  years ago (Martel et  al. 2013). This means that they are still hot (500– 600 °C), if not partially melted. Surface heat flow and heat production Part of the Pontgibaud fault zone was intersected by core drilling undertaken by the BRGM (French Geological Survey) in the 1970s. The purpose of these drillings was to characterize the mining potential of these structures. As part of this study, one of these boreholes (Peyrouses 1, location in Fig.  2) remains accessible and has been used in this study. In the past, temperature measurements have been carried out between 55 and 120 m. From these measurements, a geothermal gradient of 41 °C/km has been inferred −2 and a heat flux of 110  mW  m was estimated (IHFC database; http://www.datap ages. c om/g i s -ma p- publi s hing - pr o g r am/g i s -op en-file s /g lob a l-f rame wor k/g lob a l-he a t -flow - datab ase). This relatively high heat flow value can be explained by an elevated mantle heat flow and/or a large crustal heat production. The analysis of 660 rock samples taken from granitic units in the FMC (Lucazeau 1981) revealed an average heat production of −3 4.06 μW m . Considering an average granite thickness of 7 km (Vigneresse et al. 1999), the heat production of this granitic upper crust can only explain 26% of the heat flow measured in the Pontgibaud fault zone. The analysis of 112 samples taken from meta - −3 morphic units revealed an average heat production of 2.27 μW m . If this average value is representative for the remaining 22 km of crust (considering the crustal thickness in Duwiquet et al. Geotherm Energy (2019) 7:33 Page 7 of 29 −2 this area is 29 km, Zeyen et al. 1997), the mantle heat flow would be close to 32 mW m , twice as large as for a thermally stable Archean crust. Permeability In a crustal context, permeability and connected porosity values are generally considered −18 very low to medium, with porosity ranging from 0 to 15% and permeability from 10 −13 2 to 10  m (Brace 1984; Evans et  al. 1997; Sonney and Vuataz 2009; Ingebritsen and Appold 2012; Ranjram et  al. 2015; Walter 2016; Achtziger-Zupancic et  al. 2016, 2017). However, the orders of magnitude for permeability values and derived parameters vary greatly from one site to another and need to be addressed specifically to better charac - terize a specific site or case study (Manning and Ingebritsen 1999). Therefore, new per - meability measurements on representative samples of the various lithologies from the Pontgibaud fault zone have been performed. An additional way to constrain the fault zone system is to use geophysical data, as detailed below. Magnetotelluric data Magnetotelluric (MT) campaigns were carried out in the Pontgibaud fault zone area by TLS-Geothermics and IMAGIR between 2015 and 2017 (Ars et al. 2019). From their 3D model of electrical resistivity produced in 2016 by IMAGIR, the profile corresponding to our geological cross-section was extracted (Fig.  3). The interest of this passive elec - tromagnetic method is its sensitivity to geometries and its ability to image conductive objects in three dimensions. Results show that the obtained conductive bodies in the upper crust could be related to the Pontgibaud fault zone. Mining work on the Pontgibaud fault zone has highlighted structures with a dip of 70° E (Bouladon et al. 1961), but more recent data favors dips of ~ 60°, as detailed later. This Pontgibaud fault zone Aigueperse-Saint-Sauves fault zone Fig. 3 2D resistivity profile. The dashed lines represent the hypothetical tracing of the deep structures, constrained by surface signatures (fault orientations and dips). (Ars et al. 2019) Duwiquet et al. Geotherm Energy (2019) 7:33 Page 8 of 29 steep dip is consistent with the structures that are visible on the MT image (Fig. 3). The obtained resistivity model makes it possible to identify highly conductive bodies, one of which is located between a depth of 10 and 20 km (Fig. 3). Petrological observations on trachytes suggest that the conductivity anomaly located at a depth between 10 and 20 km can reflect the presence of silicate melt (Martel et al. 2013). At a depth of 8 and 3  km, two other anomalies located in the brittle crust could be linked to the presence of magmas, clays (smectites) and/or hydrothermal fluids. What - ever their nature, these two anomalies appear to be related to the Pontgibaud fault zone. This fault zone, which could be interpreted as a listric fault zone, seems to join the Aigueperse–Saint-Sauves fault system at a depth of 10 to 15  km (Ars et  al. 2019). Isotopic analyses of CO from the Ceyssat and Pranal resurgences at δ C (Fig.  2) indi- cate a mantle origin (Verney 2005), which is in agreement with isotopical results from the whole French Massif Central (Brauer et al. 2017). Indeed, the presence of CO at the surface can be explained by the connectivity of the Pontgibaud and Aigueperse–Saint- Sauves structures to magmatic activity. It also highlights that these two fault zones serve as pathways for hydrothermal fluid flow. Methods Field observations This approach aims to characterize the intensity of the brittle deformation as well as the general geometry of the Pontgibaud fault zone. Structural data (fault orientations and dips) collected since 2014 by TLS-Geothermics have been compiled to assess the inten- sity of deformation. These observations are key to improve the geometry and scaling of the simplified numerical model. The microstructural analysis of the samples collected within this hydrothermal zone (“Peyrouses 1” borehole). Figure  2 provides an under- standing of how fluids flow within this deformed zone. Laboratory observations and measurements In order to qualitatively and quantitatively characterize the permeability and porosity parameters of the hydrothermal zone, 2D thin-section observations, 3D X-ray micro- tomography observations and laboratory measurements of permeability and connected porosity were performed. The thin-section observations were performed on the same samples that were analyzed by X-ray micro-tomography. X-ray micro-tomography is a non-destructive imaging technique that provides a three-dimensional visualization of the analyzed samples. X-rays emitted by a tungsten filament pass through the sample. Its power depends on the thickness and nature of the radiographed element. A detector, located behind the sample, measures the attenua- tion of the X-ray beam. A series of images is then generated in different greyscales. A three-dimensional reconstruction of the object is performed using this series of images. However, the density contrast that may exist between dense mineral phases (pyrite) and low-density mineral phases (feldspar) can bias the information. Beam–Hardening–Cor- rection (BHC) limits this loss of information. This imaging technique makes it possible to obtain voxel resolutions between 25 and 4  μm. Once rebuilt, the three-dimensional volume analysis and processing can be performed using the Volume Graphics VG Studio Max software (https ://www.volum egrap hics.com/en/produ cts/vgstu dio-max.html). Duwiquet et al. Geotherm Energy (2019) 7:33 Page 9 of 29 The connected porosities and permeabilities of selected samples were measured at the Ecole et Observatoire des Sciences de la Terre (EOST) in Strasbourg (France). Measure- ments were performed on the seven samples for which the porosity was visualized using X-ray micro-tomography (including fractured samples, intact samples, altered samples, and altered and fractured samples). First, the oven-dried (in a vacuum at 40 °C for at least 48 h) 20 mm diameter and 40 mm long core samples were placed in an AccuPyc II helium pyc- nometer to calculate their skeletal (connected volume). The connected porosity (the inter - connected porosity within the sample) was then calculated using the bulk volume of the samples, measured using digital calipers. Second, the gas permeability was measured using a benchtop nitrogen permeameter (Farquharson et al. 2016; Heap and Kennedy 2016; Heap 2019). All measurements were made at room temperature and under a confining pressure of 1 MPa. Permeability was measured using either the “steady-state” or the “pulse-decay” method, depending on the permeability of the sample [methods are outlined in detail in Heap et al. (2017) and Kushnir et al. (2018)]. For the steady-state measurements, volumetric flow rates were measured (using a gas flow meter) for different values of differential pressure (upstream minus the downstream pressure; measured using a pressure transducer). For the pulse-decay measurements, the pressure decay of a known volume of gas (in the upstream circuit) was monitored as a func- tion of time. Estimate of Rayleigh number The Rayleigh number expresses the ratio between the driving forces (e.g., gravity) that pro - mote fluid circulation and the resistant forces (e.g., viscous resistance) that prevent circula - tion. It is expressed in the form: KLα�Tg Ra = , (1) Dν 2 −1 where K (m ) is the medium permeability, L (m) is the thickness of the system, α (°C ) is the thermal expansion coefficient, ∆T (°C) is the temperature difference between the −2 upper and lower limits, g is the acceleration due to gravity (m s ) and ν is the kinematic 2 −1 2 −1 fluid viscosity (m  s ). D is the thermodispersion tensor (m  s ) (e.g., Magri et al. 2016) as defined by: Φf + (1 − Φ)s D = , (2) ρf Cpf −1 −1 where Φ corresponds to the porosity, Cpf (J  kg  K ) to the fluid thermal capacity, λf −1 −1 and λs, respectively, to the fluid and solid thermal conductivity (W m  K ). The critical Rayleigh number is the value above which heat is evacuated through convec - tive circulation. For 3D vertical faults with a permeable top, and in the case of a constant fluid viscosity (not temperature-dependent), Malkovsky and Pek (2004) provided an analyt - ical relationship to calculate the critical Rayleigh number (Rcf ). This relationship depends only on the fault aspect ratio Δ: 0.8584 1.165 6.428 1.165 Rcf = + (27.1) , (3) � Duwiquet et al. Geotherm Energy (2019) 7:33 Page 10 of 29 where , where d is the thickness of the fault zone (m) and H its depth (m). Other 2H numerical and analytical simulations have shown that, for a temperature-dependent fluid viscosity, the critical permeability above which convection begins is lower than the case where the fluid viscosity is considered constant (Malkovsky and Magri 2016). According to the calculations of Malkovsky and Magri (2016), the critical fault perme- ability Kfc from which convection is likely to be initiated in the faults is calculated using Rcf /4. Numerical approach In order to quantify the effect of fault zone dip and the impact of variable permeability on fluid circulation in the fault zone, a series of numerical simulations using Comsol Multiphysics software was performed. Heat equations and Darcy’s law were coupled in a similar way to the Guillou-Frottier et al. (2013). The numerical approach was organized in two successive stages. The first stage was dedicated to the study of a theoretical fault zone and focused on the influence of fault dip and permeability on fluid flow in order to understand the processes governing deep fluid circulation. The second stage was dedi - cated to the Pontgibaud hydrothermal system, in order to verify the legitimacy of the digital crustal-scale model. Finally, the obtained results, and in particular the geothermal gradient and surface heat flow, were compared with field data, including the resistivity data from Ars et  al. (2019) (Fig.  3) which were not considered to build our numerical model. Boundary and initial conditions The average temperature at the surface was fixed at 10  °C. According to regional data (see “Geothermal setting of the Pontgibaud area” section), a constant heat flux of 35 −2 mW  m was imposed at the base of the model. The lateral borders of the model were thermally insulated. The initial conditions of the model were defined from the stationary state obtained for the pure thermal conduction regime. At t = 0, the permeabilities were assigned to the different compartments, and the fluid flow was initiated. Considering the thermal relaxation time of the magmatic chambers (see above, Martel et  al. 2013) calculations were performed up to a time of t + 15,000 years. For Darcy’s law, a no flow condition was imposed at the bottom and on both sides of the model, a surface pressure equivalent to atmospheric pressure (P = 10 Pa) was imposed at the surface and a positive pressure evolution with depth was adopted as the initial condition. Meshing and geometry Two geometrical configurations were tested. The first configuration was that of the theo - retical fault zone (Fig.  4a), and the second was the crustal-scale model applied to the Pontgibaud case study (Fig.  4b), which is constrained by observations and field meas - urements. For the smallest (i.e., for the models of the theoretical fault zone) fault zone models, the mesh consisted of 61,211 triangles with a size not larger than 5  m in the high-permeability zones and not smaller than 150 m in the low-permeability zones. The mesh was also refined at the intersections of the fault planes. For the largest model (i.e., Duwiquet et al. Geotherm Energy (2019) 7:33 Page 11 of 29 Fig. 4 Geometries studied by numerical modeling. a Theoretical fault zone geometry. b Simplified crustal‑scale geometry containing essential structural elements derived from field and geophysical data, and MT interpretation the crustal models), the mesh consisted of 52,798 triangles with a size not larger than 25 m in the high-permeability zones and not smaller than 300 m in the low-permeability zones. Variable physical properties The dynamic viscosity of the fluid is function of temperature T (°C), such that: −5 µ (T ) = 2.414 · 10 · exp (4) T + 133 according to Smith and Chapman (1983) and Rabinowicz et  al. (1998). The density of water is also a function of temperature: ρ (T ) = 1002.4 − 0.1905 · T − 0.0025 · T (5) based on a polynomial adjustment of experimental data up to 350 °C, as used by Guil- lou-Frottier et  al. (2013). For temperatures above 350  °C, the fluid density was fixed at −3 600  kg  m . This hypothesis can be considered valid since such temperatures within Duwiquet et al. Geotherm Energy (2019) 7:33 Page 12 of 29 −16 2 the fault are reached at depths where permeability is lower than 10 m , preventing buoyancy forces to be effective. The permeabilities K of the basement (m ) and K of the b f faults (m ) are the main adjustment variables. For the basement, the decrease in perme- ability with depth can be physically characterized by the following equation (Garibaldi et al. 2010): z − 800 K (z) = K exp , (6) b b where K (z) is the permeability of the depth-dependent basement, K is the permeabil- b b ity at the surface, and z is the depth below sea level. The length δ (m) characterizes the intensity of the decrease in permeability with depth, taken here at 2500 m. As the Gelles granite and the pre-Tertiary bedrock are intensely fractured, Eq. (6) was applied to these lithologies, in accordance with the models of Garibaldi et al. (2010) and Guillou-Frottier et  al. (2013). According to laboratory measurements (see below), a permeability value −16 2 of 10 m was assigned to K . In addition, to account for the uncertainty associated with upscaling the experimental values (Neuman 1994; Heap and Kennedy 2016), a large range of values has been tested in the numerical models. Within the fault zone, we considered that the permeability variation with depth varied according to Eq.  (7). However, field observations indicate that the area has a very high fracture density. In addition, the lateral variation in permeability will be very heterogeneous. Therefore, we consider that the use of Eq.  (6) is not suffi- cient. By keeping the decrease in permeability with depth, a lateral variation can be included, by using the following spatial variation: z − 800 z + x ( ) ( ) K (x, z) = K × exp × 101 + 100 × sin 2 × π × , f f (7) where K (x, z) is the space-dependent permeability of the fault and K × 201 being the f f maximum permeability value at the surface. The last term on the right side of Eq.  (7) makes the permeability alternate along a sinusoid applied to the vertical and horizontal axes of the numerical model. The term λ corresponds to the wavelength of the sinusoid. To reproduce rather fine alternations of high and low permeability, as suggested by field observations (Fig.  6), a low value of λ was chosen. Finally, the choice of constants 100 and 101 makes it possible to vary the permeability over two orders of magnitude. Numerical modeling allows us to test the effect of the unknown parameters. In order to see the effects of geometry and permeability, the physical properties on either side of the fault zone were identical and fixed (Fig. 5). During the various sim- ulations, the dip of the fault zone was tested at angles ranging from 10° to 90°. The assigned permeability to the host and fault zone varied according to Eqs. (6) and (7). −18 −16 2 The K value varied from 1 × 10 to 2 × 10 m , thus giving a maximum perme- −14 2 ability value ( K ) of 4 × 10 m , hence greater than the critical permeability esti- max −16 2 mated for the Pontgibaud fault zone. For the basement, the K value was 10 m −20 2 and decreased down to 1 × 10 m at depth. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 13 of 29 Fig. 5 Permeability imposed in the numerical models of the parametric study. The alpha angle corresponds to the dip of the fault zone (here α = 90°). The permeability of the basement and fault will decrease with depth according to Eq. 6 (Garibaldi et al. 2010). Within the fault zone, permeability will also vary laterally (Eq. 7) to match field observations Results Structure characterization With 126 new measurements of fault dips and orientations, three families of faults have been identified. The first family comprises 67° ESE dipping normal N 10° trending vein beam faults (Fig.  6a). With a surface thickness of 3 km, this normal fault zone has low- ered the topography by 300  m eastwards. In addition, over the study area, the Sioule River surprisingly flows north–south over more than 10 km, and partly borders the vein network, before becoming meandering again. This north–south direction is also consist - ent with the main direction of the Limagne fault, which is located a little further east. Due to field data and geophysical arguments (see “Magnetotelluric data ” section), the general geometry of the Pontgibaud fault zone could be regarded as a listric type. The second and third fault families, both steep, intersect the Pontgibaud fault zone. The second family is formed of mainly dextral N 30° trending faults (Fig.  6b), while the third family consists of sinistral faults that trend N 115° (Fig. 6c). To summarize all field observations and data from the literature, a compiled geologi - cal cross-section is shown in Fig.  7. Its location corresponds to the geophysical profile (Fig. 3) and it crosses the villages of Prondines, Ceyssat, and ends at the Limagne basin. This cross-section illustrates all the geological and geophysical features detailed in the previous sections. Thus, the proposed listric geometry of Pontgibaud fault zone is asso - ciated with two other vertical fault families. The presence of springs within these struc - tures is an indicator of present-day fluid circulation. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 14 of 29 Fig. 6 Observations of the different lithologies and fault families present in the study area. a′, c′ Metamorphic series. b′ Granite. Schmidt’s projection on the lower hemisphere (a–c). Isodensity obtained with Stereonet ( TLS‑ Geothermics, modified) Fig. 7 Synthetic and simplified geological cross‑section of the Pontgibaud area Laboratory observations and analyses Qualitative approach: thin‑section observations Thin-section observations show fracture-rich and highly altered facies (Fig.  8). These fractures are mainly filled with quartz, feldspar, phyllosilicates, ferriferous, pyrite and oxide minerals. If at the heart of these fractures the circulation of fluids appears diffi - cult, the walls seem to be favorable for the circulation of fluids (Fig.  8c). Fluid circulation within the fractures promoted the crystallization of secondary mineral phases. Turma- line (Fig.  8d), is observed in large quantities in all thin-sections. Fluid circulation also caused significant mineral dissolution (Fig.  8b). This dissolution can be accompanied by a loss of volume, thus promoting the formation of pores within the matrix (Launay 2018). These initial observations show that fractures related to tectonic history served as pathways for (paleo)fluids and produced haloes of alteration. These haloes were subject to significant mineral dissolution resulting in the formation of micro-porosity. Fractures and micro-porosity represent the available space for modern-day fluid circulation. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 15 of 29 Fig. 8 Microscopic observations of thin‑section from the hydrothermalized zone. a Anatexite rich in veins. b Fractured anatexite. c Mineralized breach. d Porous breach. (Po) Porosity. (Qz) Quartz. ( Void) Void. (Fr) Fracture. ( Tur) Turmaline. (Ox) Oxide. (Sr) Sericites. (Kfs) Potassic feldspar However, the ability of fluids to flow through the hydrothermal system will depend on the connectivity of these fractures and micro-porosity. Three-dimensional imaging is required to qualitatively estimate the connectivity of these microstructural elements. Qualitative approach: X‑ray micro‑tomography The results of the reconstruction and image processing are given in Fig.  9. Figure  9a shows the full sample with a voxel resolution of 25  μm. It is characterized by a matrix (grey in Fig.  9a) and by fractures (yellow and red in Fig.  9a). Two phases are present within these fractures. The red phase characterizes the mineralization of the fracture. Thin-section observations on this sample showed that this fracture is filled with lead sulphides. The yellow phase represents the planar porosity (i.e., fractures). The pro - cessing of this image (Fig.  9a) made it possible to isolate these fractures (Fig.  9b) and qualitatively characterize the connectivity of the planar porosity. The three-dimensional visualization (Fig.  9b) shows a high density of fractures. These fractures are distributed heterogeneously in space. This three-dimensional arrangement promotes connectivity of the pores located at the fracture walls (yellow in Fig.  9b). Because a visualization of the matrix porosity (Fig.  9b) was not possible at this resolution, due to the characteristic pore size, a smaller volume of this sample was analyzed with a resolution of 4 μm in an attempt to detect and segment the potentially connected pores. The result of the analysis and processing of the sample is given in Fig. 9c. Coupling thin-section (2D) and micro- tomography (3D) observations highlight altered planar (i.e., fractures) and matrix (i.e., pores) porosities partially filled with secondary minerals. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 16 of 29 Fig. 9 Three‑ dimensional observations of the planar and matrix porosity of one sample taken from the Pontgibaud fault zone. The red color of figures a, b corresponds to mineralization (lead sulphides). The yellow color of figures a, b corresponds to the planar porosity located on the fracture walls. The yellow color of figure c could correspond to the matrix porosity. The resolution of figures a–c is, respectively, 25, 25 and 4 μm Because this small-scale porous network could play a crucial role in determining the fluid circulation within the Pontgibaud hydrothermal zone, we performed direct meas - urements of connected porosity and permeability in the laboratory. Quantitative approach: porosity and permeability measurements −18 −13 2 The measured permeability values varied between 2 × 10 and 7 × 10 m while the connected porosity values varied between 6 and 22% (Fig. 10). In general, permeability increases as connected porosity increases (Fig. 10). The high - est porosity and permeability were measured for the fractured and altered sample, and the lowest porosity and permeability were measured for the intact and unaltered sample (Fig. 10). To summarize, the three above-described results suggest that permeability is increased by the presence of system-lengthscale fractures and alteration (Fig.  10). The relatively high permeability of the measured samples (the permeability of granite can be as low −22 2 as 10 m ; Meredith et  al. 2012), suggest that fluid flow within the matrix could also be significant, alongside large-scale fractures and fracture networks. Hence, after hav - ing determined the first unknown parameter in “Structure characterization ” section Duwiquet et al. Geotherm Energy (2019) 7:33 Page 17 of 29 Fig. 10 Connected porosity and permeability measurements on sample taken from the Peyrouse 1 borehole (identical to those used for qualitative observations). The non‑fractured and unaltered samples have the lowest permeability and connected porosity value. Conversely, the altered and fractured sample has the highest value of permeability and connected porosity. These purely qualitative degrees of alteration and fracturation are based on macroscopic, microscopic observations and three‑ dimensional visualization by X‑ray micro ‑tomography analysis (fault geometry), we identify here the second one (permeability variation). To better understand how fault geometry and permeability variation influence hydrothermal fluid flow and thus control the establishment of thermal anomalies, a numerical approach is required. Numerical approach Rayleigh number analysis Using the Pontgibaud fault zone parameters (H = 15,000 m and h = 3000 m) and Eq. (3), the critical Rayleigh number Rcf /4 equals 21.3. The measured permeability varies from −18 −13 2 10 to 10 m (Fig.  10), implying a variation of the Rayleigh number (Eq.  1) from 0.018 to 1800. We can thus deduce the critical fault permeability above which thermal convection can occur: Rcf · D · ν −15 2 (8) K = = 1.2 × 10 m , fc g · α · �T · L −4 −1 ◦ −6 2 −1 −6 2 −1 where α = 2 × 10 K , T = 450 C , D = 0.76 × 10 m s , and ν = 10 m s . u Th s, considering the selected parameters, it seems reasonable to consider a convective thermal regime within the Pontgibaud fault zone, at least in the case when the perme- −15 2 ability is greater than 1.2 × 10 m . This threshold value must, however, be considered with caution since the fluid properties considered here correspond to those for room temperature conditions. Numerical simulations of a theoretical fault zone More than 200 simulations were performed, with the numerical set-up of Fig.  5. The results are summarized in Fig. 11. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 18 of 29 Fig. 11 Graphical representation of all numerical models of temperature anomalies on the theoretical fault zone performed (n = 216). The values of positive thermal anomalies are shown in red. Their depths are represented in brown The images in Fig.  11 represent temperature anomalies (∆T), corresponding to abnormally hot and abnormally cold temperatures compared to an undisturbed (con- ductive) thermal regime. The amplitude of the temperature anomalies is different for each case, but numerical values are indicated for some of them. In order to account for both the permeability of the fault and of the basement, the R ratio is defined by: max R = . (9) This ratio compares the maximum permeability of the fault and of the basement −16 2 (the permeability K being set at 10 m ). Other K values have been tested, but b b 0 0 higher values are not realistic and lower values do not change the results detailed below. As it can be seen in Fig. 11, for R = 200, as the dip of the fault zone increases from 10° to 90°, the temperature values increase from 12 to 85.4 °C (red numbers in Fig. 11). Thus, for a fixed R ratio, the highest thermal anomaly value will be for a vertical fault zone. The depth of the thermal anomaly also varies as a function of fault dip. By increasing the fault dip from 10° to 90°, the depth of the thermal anomaly (brown numbers in Fig. 11) decreases from 3.6 to 0.7 km. Thus, for a fixed R value, vertical structures will result in the largest thermal anomalies at the shallowest depths. By observing the variations in the temperature anomalies, we can see that when R increases for a fixed angle of 60°, the values of the thermal anomalies are 68.6  °C, 82.5 and 77.4  °C. Indeed, an increase in the R ratio does not show a significant increase in the temperature anomaly. However, the depth of these same anomalies decreases from 2.7 to 0.5 km. So for a fixed dip, when R ratio is high, the thermal anomaly will be at a maximum at a shallower depth. Therefore, a fault zone with a large dip and a high R ratio will result in the largest thermal anomalies at the shallowest depths. While this repre- sentation characterizes the distribution and intensity of positive thermal anomalies, it does not characterize the processes at their origin. Based on the results of these same numerical simulations, Fig.  12 graphically summarizes the morphology of the numeri- cally obtained isotherms as a function of fault dip angle, value of positive thermal anom- aly and the R ratio. Three main zones can be described: Duwiquet et al. Geotherm Energy (2019) 7:33 Page 19 of 29 Fig. 12 Diagram of convective regimes. Grey lines represent isotherm morphology. The dashed black line separate three area. In blue, the area of unicellular weak type convection zone. In orange, the area of unicellular medium‑type convection zone. In red, the area of bicellular strong‑type convection zone. The dashed grey contour on the right of the diagram illustrates the possible range for the Pontgibaud fault zone (i) For R values between 1.9 and 250 (blue area in Fig.  12) the morphology of iso- therms is slightly deformed. This low deformation occurs at dip values below 50°. The value of the resultant thermal anomaly will be lower than 65 °C. These results are consistent with Fig.  11 since the lower the dip of the structure and the lower the R ratio, the lower the thermal anomaly produced. The black arrows show that the circulation of fluids is controlled by a single long-wavelength convection cell. This domain will be called the “unicellular weak-type convection zone”. (ii) For R values between 40 and 400 with a dip value between 10 and 90° (orange area in Fig.  12), the isotherms show moderate deformation. The range of values of the thermal anomalies produced is between 10 and 120 °C. In the same way, the arrows show that fluid circulation is controlled by moderate-wavelength convection cells. This domain will be called the “unicellular medium-type convection zone”. (iii) Finally, for R values between 120 and 400 and dips greater than 50° (red area in Fig. 12) the isotherms are significantly deformed. The value of the thermal anomaly is higher than 67  °C. These results also show consistency with Fig.  11 since the higher the dip and R ratio, the greater the thermal anomaly at shallow depth pro- duced. The black arrows show that the circulation of fluids is controlled by two small-wavelength convection cells. This domain will be called the “bicellular strong-type convection zone”. The results of this parametric study should be applicable to hydrothermal systems in a crustal context. To sum up, for unicellular convection, the values of the thermal anomalies produced are between 10 and 120 °C. For bicellular convection, the values of the thermal anomalies produced are between 70 and 120  °C. In the next section, we apply these results to the case study of the hydrothermal system at Pontgibaud. According to Fig.  12, the Pontgibaud fault zone should fall somewhere within these two domains. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 20 of 29 Numerical simulations of the Pontgibaud fault zone The geometry of the numerical model for Pontgibaud (Fig.  13) is directly based on the geological model (Fig.  5) and is therefore constrained by a detailed set of geological data and observations. The physical properties, integrated into the numerical models, are summarized in Appendix A. In accordance with the thermal relaxation time of the magmatic chambers calculated on the basis of a thermo-barometric study (Martel et al. 2013), the results of our numerical simulations are shown up to time t + 15,000 years. Figure  13 presents the results of the numerical simulations with the temperature field (colors) and isotherms in white. The stable temperature pattern shown in Fig.  13a −15 2 is obtained for K value up to 1.4 × 10 m . For these K values, the isotherms f f max max are not deformed, but are disturbed around the cooling magma chambers. At time Fig. 13 Maps of temperatures (a), permeabilities (b), isotherms (c), and thermal anomalies (d) of the Pontgibaud region (at time t + 15,000 years). In accordance with Fig. 11, color of boxes contours correspond to one convective regime: blues boxes correspond to areas of unicellular weak‑type convection zone, oranges boxes correspond to unicellular medium‑type convection zone, and reds boxes correspond to bicellular strong‑type convection zone Duwiquet et al. Geotherm Energy (2019) 7:33 Page 21 of 29 t + 15,000 years, the cooling of the magma chambers does not disturb any isotherm in the Pontgibaud fault zone. Figure  13b represents the permeability field centered on the red box of Fig. 13a. The deformed isotherms shown in Fig.  13c, d are located at the inter- section between the Pontgibaud listric geometry and the vertical faults (see black box, −15 −14 2 Fig.  13a). For K values between 1.6 × 10 and 1.6 × 10 m , the deformation of max the isotherms can be observed (Fig.  13c). In the terms previously used, the morphol- ogy of the isotherms changes from areas of weak convection of unicellular type (blue −15 2 box, K = 1.6 × 10 m ) to areas of strong convection of bicellular type (red box, max −14 2 K = 1.6 × 10  m ). This representation therefore shows that, when permeability is max increased within the Pontgibaud fault zone, isotherms show an increasing disturbance. Figure  13d shows the amplitudes and locations of the obtained thermal anomalies. These anomalies are represented as a function of R ratio, the value of the maximum tem - perature anomaly (∆T ) and its depth. As it can be seen in Fig. 13d, positive and nega- max tive thermal anomalies are present in the Pontgibaud system. Positive thermal anomalies are located at the top of the fault zone, while negative thermal anomalies are located at the base of the fault zone. The depth of the positive thermal anomaly decreases for an increasing R ratio (see dashed line in Fig.  13d). This latter observation is consistent with the parametric study (Figs.  11, 12). For R ratios between 42 and 160, the value of the thermal anomaly does not show significant variation. Overall, these anomalies are between 57 and 66 °C. On the other hand, when R = 16, the value of the thermal anom- aly is only 24 °C. This lower value can be explained by a low deformation of the isotherms. At R = 42 and above, the isotherms show an increasing deformation up to R = 160. Considering isotherms with a relatively large deformation, the value of the positive thermal anom- aly no longer varies significantly. Thus, the results of numerical simulations carried out within the Pontgibaud hydrothermal system show that the deformation of the isotherms −15 2 is significant from a K value of 1.6 × 10 m (in accordance with the previous Ray- max leigh number analysis) and a R ratio of 16. Combining field data, laboratory measurements and hydrothermal models Finding the right R and K combination for Pontgibaud fault zone The comparison of the results obtained by modeling with the data measured in the field makes it possible to assess the consistency of the obtained results. Within the Pontgibaud −2 fault zone, the surface heat flux is 110 mW m and the geothermal gradient is 41 °C/km (Vasseur et al. 1991, International Heat Flow Commission). This comparison will make it possible to estimate the permeability value for which the field measurements and simu - lation results are consistent. Among the numerous tested models, it turns out that only some combinations precisely reproduced the above-described thermal features. The −14 2 result for the simulation with a K = 1.6 × 10 m , show a surface heat flux of 115 max −2 −1 −14 2 mW m and a geothermal gradient of 39 °C km . Thus, for K = 1.6 × 10 m , the max results of the numerical models reproduce the heat flux and geothermal gradient values measured near the surface. The numerical models show that at a time t + 10,000  years −2 the heat flux and the geothermal gradient are 105 mW m and 33 °C/km, respectively. −2 At t + 20,000 years we obtained 127 mW m and 41 °C/km, respectively. 0 Duwiquet et al. Geotherm Energy (2019) 7:33 Page 22 of 29 Combining geophysical data and numerical results Figure  14a, b and c compare the temperature anomaly map of the Pontgibaud region with the measured resistivity profile extracted at the exact coordinates of the geological cross-section. The map (Fig.  14a) shows positive thermal anomalies within the Pontgibaud fault zone, at the magma chambers and at the base of the sedimentary infill of the Limagne basin (right part of the model). A negative thermal anomaly is also present at the base of the Pontgibaud fault zone. The MT image shows various anomalies of low resistivity at 3, 8 and 15 km depth (Fig.  14b). By superimposing the two images (Fig.  14a, b), we can see that the temperature anomalies located at 3, 8 and 15 km depth are superimposed on the low-resistivity anomalies located at 3, 8 and 15 km depth. At 15 km depth, the thermobarometric study by Martel et al. (2013) shows the pres- ence of a magma chamber. At 3 and 8  km depth, Fig.  14c shows that anomalies of low resistivity are located in areas where fluids could circulate. Thus, the numerical simula - −14 2 tion imposing a permeability K = 1.6 × 10 m at the Pontgibaud fault zone shows max a striking consistency with the resistivity profile, even if the link between resistivity and Fig. 14 a Thermal anomalies maps of Pontgibaud area. The black arrows represent the direction of fluid flow. b Resistivity model previously described (Fig. 3). c Comparison of the numerical simulation result affecting −14 −2 a K permeability of 1.6 × 10 m at the Pontgibaud fault zone, with the resistivity profile. d Fluid flow max velocity within the Pontgibaud fault zone Duwiquet et al. Geotherm Energy (2019) 7:33 Page 23 of 29 Fig. 15 Numerical simulation results for a maximum permeability imposed on the Pontgibaud fault zone of −14 −2 1.6 × 10 m temperature may not be so obvious. Nevertheless, isotherms deformation is correlated with high fluid flow velocity (Fig. 14d). Discussion Multi‑disciplinary approach This study aims to understand how crustal-scale fault zones can generate a viable geo - thermal resource. The overall understanding of an outcropping fault system inte - grates many parameters that can be studied using a multi-disciplinary approach. Thus, structural observations and measurements, 2D and 3D thin-section and X-ray micro- tomography observations, as well as connected porosity and permeability laboratory measurements, have highlighted two essential variables: the fault dip angle and perme- ability of a fault zone. Numerical modeling has made it possible to test the effect of these two variables on the intensity and depth of positive temperature anomalies, but also to characterize the deformation of the resultant isotherms. The comparison of the results of the numerical simulations at crustal scale with the field data shows a permeability −14 2 of the Pontgibaud fault zone K = 1.6 × 10 m . For this value, the ratio R = 160 max (Figs.  11, 12). The value of the thermal anomaly varies between 68 and 70  °C and has depths between 2.7 and 1.7  km (Figs.  11, 12). Our results show a thermal anomaly of 66 °C within the Pontgibaud fault zone at a depth of 1.8 km (Fig.  15). It is worth men- −14 2 tioning that the permeability value of 1.6 × 10 m is in accordance with previous Duwiquet et al. Geotherm Energy (2019) 7:33 Page 24 of 29 studies where temperature measurements and surface heat flow values were numerically reproduced with similar permeability values (e.g.,  Person et  al. 2012; Guillou-Frottier et al. 2013; Roche et al. 2018a, b). Origin of porosity Thin-sections observations (Fig.  8) revealed open fractures and fractures completely sealed by mineral precipitation. A fault can thus behave as a pathway for, or as a barrier to, fluid circulation (Bense et  al. 2013). Quartz, feldspar and alteration products such as sericites are well known to form around hydrothermal sites and near faults (Bruhn et  al. 1994; Sausse 1998; Caine et  al. 2010). Filling fractures with cement or alteration products can reduce their permeability (e.g., Griffiths et al. 2016; Mordensky et al. 2018). When fractures transmit fluids, their immediate surroundings can be subject to signifi - cant mineral dissolution (potassic feldspar, biotite, plagioclases), forming matrix poros- ity around these fractures. If the minerals likely to be altered are initially connected in space, then their alterations will produce connected pores that can promote the circula- tion of fluids (e.g., Farquharson et al. 2018). Direct connected porosity measurements on our samples with no apparent fractures (at the lengthscale of the 40  mm-long sample) showed values between 6 and 22%. Indeed, the permeability of the altered samples was higher than that of the unaltered samples (Fig.  10). The porosity values measured here are consistent with the study by Sardini et al. (1997) who quantified the connectivity of the minerals likely to be altered. Mercury porosity between 6 and 20% was measured on samples of altered, non-fractured granite. Considering the capacity of fluids to circulate through fractures but also within the matrix itself, it is thus possible to consider Pon- tgibaud hydrothermal system as a reservoir with double matrix-fracture permeability. This kind of reservoir has already been observed in the altered granite of the Soultz- Sous-Forêts geothermal reservoir (Genter 1989; Ledésert 1993). Controlling fluid flow factors within the Pontgibaud hydrothermal system Although we did not consider the present-day stress field, the value of thermal anom - alies and their depths are the result of convective heat transfer processes. López and Smith (1995) showed that this thermal regime is directly related to the permeability con- trast between the fault and its host rocks. In addition to the results of López and Smith (1995), in which the R ratio was studied, we have explored here the effect of the dip angle of the fault zone on fluid circulation. Putative structures sufficiently permeable to allow fluid circulation are commonly described for hydrothermal systems (Forster and Evans 1991; Bruhn et al. 1994; Belgrano et al. 2016) and are used in many numerical models (Forster and Smith 1989; Mckenna and Blackwell 2004; Garibaldi et al. 2010; Guillou-Frottier et al. 2013; Magri et al. 2016; Taillefer et al. 2018a, b). The results of the parametric study presented herein, and their comparison with the Pontgibaud natural system, show that when the fault zone is verti- cal, the value of thermal anomaly will be high and at shallow depth. This depth is also controlled by the ratio between the fault and its host permeability: when the R ratio is high, thermal anomaly will be at a shallower depth. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 25 of 29 Permeability measurements were carried out on samples from the Pontgibaud fault −18 −13 2 zone (Peyrouses 1 drilling). The results show values between 2 × 10 and 7 × 10 m . The permeabilities of the rocks measured in this study are relatively high compared to laboratory measurements on intact granite samples (e.g., Brace et al. 1968; Kranz et al. 1979; Meredith et al. 2012). The maximum permeability value for the basement assigned −16 2 to the ratio R has been set at 10 m . This value is compatible with the permeabil - ity values extracted from the laws of permeability decrease proposed by Ingebritsen and Appold (2010, 2012). Conclusions The aim of this study was to better understand the organization and controlling factors of a hydrothermal system associated with crustal faults, with an application to the Pon- tgibaud hydrothermal system. Based on results from field observations, microscopic observations and laboratory analyses, the Pontgibaud fault zone has the characteristics of a dual matrix-fracture permeability reservoir. The results also highlighted parameters that are essential for understanding a hydrothermal system: the dip angle of the fault and its permeability. Numerical 2D simulations, combining fluid flow and the heat equation, have tested the effect of these parameters on the value and distribution of positive ther - mal anomalies within the fault zone. Based on more than 200 numerical simulations, the results of this parametric study showed that (i) a vertical structure provides the largest thermal anomaly at the shallowest depth; (ii) the depth of positive thermal anomaly is shallow for high R ratio values. Although the model is quite simple, it well captures the convective behavior. These generic results were then compared to a natural case study, the hydrothermal system at Pontgibaud. Numerical simulations performed at a crustal scale, and integrat- ing observations and field measurements, showed similar results to those of the afore - mentioned parametric study. However, a 3D study should bring new insights into the Pontgibaud hydrothermal system and its associated thermal anomalies. The results of fluid circulation simulations within the Pontgibaud fault zone showed −14 2 that a maximum fault permeability value of 1.6 × 10 m allows reproduction of the surface data. Considering the economic potential, a temperature of 150 °C at a depth of 2500  m represents an interesting target for high-temperature geothermal exploitation. Finally, a new geothermal play, crustal fault zones was identified and intensively explored and analyzed. A methodology has been developed to constrain this specific hydrother - mal system, and to target economically interesting areas. Furthermore drilling planned by operators should make it possible to verify and validate all or parts of this approach. Abbreviations BRGM: Bureau de Recherches Géologiques et Minières; CFZ: Crustal fault zone; EOST: Ecole et Observatoire des Sciences de la Terre; FMC: French Massif Central; ISTO: Institut des Sciences de la Terre d’Orléans; IMAGIR: https ://www.imagi r.eu/; TLS: TLS‑ Geothermics (https ://www.tls‑geoth ermic s.fr). Acknowledgements The financial support was provided by TLS‑ Geothermics and supported by ISTO and BRGM. Funding for open access was generously offered by the Helmholtz Centre for Environmental Research‑UFZ, Helmholtz Centre Potsdam‑ GFZ German Research Centre for Geosciences, the Karlsruhe Institute of Technology (KIT ), and LABEX Grant ANR‑11‑LABX ‑0050_G‑ EAU‑ THERMIE‑PROFONDE (this research therefore benefited from state funding managed by the Agence National de la Recherche (ANR) as part of the “Investissements d’avenir” program)”. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 26 of 29 Authors’ contributions This work corresponds to M.Sc. thesis of HD. LA analyzed, interpreted and retrieved the field data. LGF contributed to numerical modeling and participated in writing the manuscript. MJH helped HD to perform the permeability measure‑ ments. MB provided a very large amount of data and participated in the various field campaigns. All authors read and approved the final manuscript. Funding The financial support was provided by TLS‑ Geothermics. Availability of data and materials All data are available in Figs. 1, 2, 4, 8 and Appendix A. Competing interests The authors declare that they have no competing interests. Author details 1 2 3 Université d’Orléans, ISTO, UMR 7327, 45071 Orléans, France. BRGM, 45060 Orléans, France. ISTO, UMR7327, Université d’Orléans, CNRS, BRGM, 45071 Orléans, France. TLS‑ Geothermics, 92 Chemin de Gabardie, 31200 Toulouse, France. Géophysique Expérimentale, Institut de Physique de Globe de Strasbourg, UMR 7516, CNRS, Université de Strasbourg/ EOST, 5 Rue René Descartes, 67084 Strasbourg Cedex, France. Appendix Appendix A Rock properties Symbols Magmatic Calcareous‑ Metamorphic Gelles Carboniferous Migmatitic– Mantle Unit chamber clay unit granite plutonic rocks gneissic complex complex Porosity Ф 1.0 10 5 5 5 3 2 % Perme‑ k Variable Variable Variable Variable Variable Variable Vari‑ m ability able Thermal λ 1.3 2.2 2.5 3.25 2.8 2.5 3 W/m K con‑ duc‑ tivity Heat Cp 800 880 800 800 800 800 1000 J/(kg K) capac‑ ity −3 Heat A 0.5 1.7 2.1 3.5 3.5 2 0.5 µW m pro‑ duc‑ tion Bulk ρ 3300 2650 2700 2600 2600 3000 3300 kg/m den‑ sity Fluid properties Symbols Values Unit Bulk density ρ Variables kg/m Thermal conductivity λ 0.6 W/m K Heat capacity Cp 4200 J/(kg K) Dynamic viscosity µ Variables Pa s Properties of these parameters have been recovered in the following works: Mareschal and Jaupart (2011), Lucazeau et al. (1984), Hasterok et al. (2017), Taillefer (2017), Guil- lou-Frottier et  al. (2013), Launay (2018), McKenna and Blackwell (2004), Smith and Chapman (1998), Rabinowicz et al. (1998). Duwiquet et al. Geotherm Energy (2019) 7:33 Page 27 of 29 Received: 21 May 2019 Accepted: 24 October 2019 References Achtziger‑Zupančič P, Loew S, Hiller A, Mariethoz G. 3D fluid flow in fault zones of crystalline basement rocks (Poehla‑ Tellerhaeuser Ore Field, Ore Mountains, Germany). Geofluids. 2016;16(4):688–710. Achtziger‑Zupančič P, Loew S, Mariéthoz G. A new global database to improve predictions of permeability distribution in crystalline rocks at site scale. J Geophys Res. 2017;122(5):3513–39. Ars JM, Tarits P, Hautot S, Bellanger M, Coutant O, Maia M. Joint inversion of gravity and surface wave data constrained by magnetotelluric: application to deep geothermal exploration of crustal fault zone in felsic basement. Geothermics. 2019;80:56–68. Barton C‑A, Zoback M ‑D, Moos D. Fluid flow along potentially active faults in crystalline rock. Geology. 1995;23(8):683–6. Baujard C, Genter A, Dalmais E, Maurer V, Hehn R, Rosillette R, Schmittbuhl J. Hydrothermal characterization of wells GRT‑1 and GRT ‑2 in Rittershoffen, France: implications on the understanding of natural flow systems in the rhine graben. Geothermics. 2017;65:255–68. Belgrano T‑M, Herwegh M, Berger A. Inherited structural controls on fault geometry, architecture and hydrothermal activ‑ ity: an example from grimsel pass, Switzerland. Swiss J Geosci. 2016;109(3):345–64. Bellanger M. High temperature geothermal resources of crustal fault zones: a dedicated approach. In: 79th EAGE confer‑ ence and exhibition, June 2017. https ://doi.org/10.3997/2214‑4609.20170 1771. Bense V, Person M, Chaudhary K, You E, Cremer Y, Simon N‑S. Thermal anomalies indicate preferential flow along faults in unconsolidated sedimentary aquifers. Geophys Res Lett. 2008;35(24):1025. Bense V‑F, Person M ‑A. Faults as conduit ‑barrier systems to fluid flow in siliciclastic sedimentary aquifers. Water Resour Res. 2006;42:5. Bense V, Gleeson T, Loveless S, Bour O, Scibek J. Fault zone hydrogeology. Earth Sci Rev. 2013;127:171–92. Bertrand L. Etude des réservoirs géothermiques développés dans le socle et à l’interface avec les formations sédimen‑ taires. Dissertation, Université de Lorraine. 2016. Bouladon JP, Picot J‑ J, Sainfeld P. Le faisceau filonien de Pontgibaud. Bull BRGM. 1961;1:1–41. Brace A‑ W. Permeability of crystalline rocks: new in situ measurements. J Geophys Res. 1984;89(B6):4327–30. Brace W, Walsh J‑B, Frangos W ‑ T. Permeability of granite under high pressure. J Geophys Res. 1968;73(6):2225–36. Bril H, Bonhomme M‑ G, Marcoux E, Baubron J‑ C. Ages K/Ar des minéralisations de Brioude‑Massiac ( W ‑Au‑As‑Sb; Pb ‑Zn), Pontgibaud (Pb‑Ag; Sn), et Labessette (As‑Pb ‑Sb ‑Au): place de ces districts dans l’évolution géotectonique du Massif central français. Miner Deposita. 1991;26(3):189–98. Brockamp O, Schlegel A, Wemmer K. Complex hydrothermal alteration and illite K‑Ar ages in Upper Visean molasse sedi‑ ments and magmatic rocks of the Variscan Badenweiler‑Lenzkirch suture zone, Black Forest, Germany. Int J Earth Sci. 2015;104(3):683–702. Brogi E, Fulignati P. Tectonic control on hydrothermal circulation and fluid evolution in the pietratonda‑poggio peloso (southern Tuscany, Italy) carbonate‑hosted sb ‑mineralization. Ore Geol Rev. 2012;44:158–71. Bruhn R‑L, Parry W ‑ T, Yonkee W‑A, Thompson T. Fracturing and hydrothermal alteration in normal fault zones. Pure Appl Geophys. 1994;142(3–4):609–44. Caine J‑S, Evans J‑P, Forster C‑B. Fault zone architecture and permeability structure. Geology. 1996;24(11):1025–8. Caine J‑S, Bruhn R‑L, Forster C‑B. Internal structure, fault rocks, and inferences regarding deformation, fluid flow, and mineralization in the seismogenic stillwater normal fault, dixie valley, nevada. J Struct Geol. 2010;32(11):1576–89. Edel J‑B, Maurer V, Dalmais E, Genter A, Richard A, Letourneau O, Hehn R. Structure and nature of the Palaeozoic base ‑ ment based on magnetic, gravimetric and seismic investigations in the central Upper Rhinegraben. Geotherm Energy. 2018;6(1):13. Evans J‑P, Forster C‑B, Goddard J‑ V. Permeability of fault‑related rocks, and implications for hydraulic structure of fault zones. J Struct Geol. 1997;19(11):1393–404. Fairley J‑P. Modeling fluid flow in a heterogeneous, fault ‑ controlled hydrothermal system. Geofluids. 2009;9(2):153–66. Famin V, Philippot P, Jolivet L, Agard P. Evolution of hydrothermal regime along a crustal shear zone, Tinos Island, Greece. Tectonics. 2004. https ://doi.org/10.1029/2003T C0015 09. Farquharson J‑I, Wild B, Kushnir A‑R, Heap M ‑ J, Baud P, Kennedy B. Acid‑induced dissolution of andesite: evolution of permeability and strength. J Geophys Res. 2018. https ://doi.org/10.1029/2018J B0161 30. Farquharson J‑I, Heap M ‑ J, Lavallée Y, Varley N‑R, Baud P. Evidence for the development of permeability anisotropy in lava domes and volcanic conduits. J Volcanol Geoth Res. 2016;323:163–85. Faulds J, Coolbaugh M, Bouchot V, Moek I, Oguz K. Characterizing structural controls of geothermal reservoirs in the Great Basin, USA, and western turkey: Developing successful exploration strategies in extended terranes. In: Pro‑ ceedings world geothermal congress Bali, Indonesia, April 25–30. 2010. Faulkner D, Jackson C, Lunn R, Schlische R, Shipton Z, Wibberley C, Withjack M. A review of recent developments con‑ cerning the structure, mechanics and fluid flow properties of fault zones. J Struct Geol. 2010;32(11):1557–75. Faure M, Grolier J, Pons J. Extensional ductile tectonics of the Sioule metamorphic series ( Variscan French Massif Central). Geol Rundsch. 1993;82(3):461–74. Fernandez A. La série cristallophyllienne et les granites de la région de Pontgibaud (Puy de Dôme) Massif Central Fran‑ çais. Dissertation, Université de Clermont Ferrand. 1969. Forster C, Smith L. The influence of groundwater flow on thermal regimes in mountainous terrain: a model study. J Geophys Res. 1989;94(B7):9439–51. Forster C‑B, Evans J‑P. Hydrogeology of thrust faults and crystalline thrust sheets: results of combined field and modeling studies. Geophys Res Lett. 1991;18(5):979–82. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 28 of 29 Garibaldi C, Guillou‑Frottier L, Lardeaux J‑M, Bonté D, Lopez S, Bouchot V, Ledru P. Thermal anomalies and geologi‑ cal structures in the Provence basin: implications for hydrothermal circulations at depth. Bull Soc Géol France. 2010;181(4):363–76. Genter A. Géothermie Roches Chaudes Sèches: le granite de Soultz‑sous‑Forêts. (Bas‑Rhin, France). Dissertation, Univer ‑ sité d’Orléans. 1989. Granet M, Stoll G, Dorel J, Achauer U, Poupinet G, Fuchs K. The Massif Central (France). New constraints on the geody‑ namical evolution from teleseismic tomography. Geophys J Int. 1995a;121:33–48. Granet M, Wilson M, Achauer U. Imaging a mantle plume beneath the French Massif Central. Earth Planet Sci Lett. 1995b;136:281–96. Grasby S‑E, Hutcheon I. Controls on the distribution of thermal springs in the southern cordillera. Can J Earth Sci. 2001;38(3):427–40. Griffiths L, Heap M ‑ J, Wang F, Daval D, Gilg H‑A, Baud P, Genter A. Geothermal implications for fracture ‑filling hydrother ‑ mal precipitation. Geothermics. 2016;64:235–45. Gudmundsson A. Active fault zones and groundwater flow. Geophys Res Lett. 2000;27(18):2993–6. Guillou‑Frottier L, Carre C, Bourgine B, Bouchot V, Genter A. Structure of hydrothermal convection in the upper Rhine graben as inferred from corrected temperature data and basin‑scale numerical models. J Volcanol Geoth Res. 2013;256:29–49. Hasterok D, Webb J. On the radiogenic heat production of igneous rocks. Geosci Front.2017;8(5):919–40. Heap M‑ J. The influence of sample geometry on the permeability of a porous sandstone. Geosci Instrum Methods Data Syst. 2019;8(1):55–61. Heap M‑ J, Kennedy B‑M. Exploring the scale ‑ dependent permeability of fractured andesite. Earth Planet Sci Lett. 2016;447:139–50. Heap M‑ J, Violay M, Wadsworth F‑B, Vasseur G. From rock to magma and back again: the evolution of temperature and deformation mechanism in conduit margin zones. Earth Planet Sci Lett. 2017;463:92–100. Haines S, Lynch E, Mulch A, Valley J‑ W, van der Pluijm B. Meteoric fluid infiltration in crustal‑scale normal fault systems as indicated by δ18O and δ2H geochemistry and 40Ar/39Ar dating of neoformed clays in brittle fault rocks. Litho‑ sphere. 2016;8(6):587–600. Ingebritsen S‑E, Manning C‑E. Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism. Geofluids. 2010;10(1–2):193–205. Ingebritsen S‑E, Appold M ‑S. The physical hydrogeology of ore deposits. Econ Geol. 2012;107(4):559–84. Jolivet L, Famin V, Mehl C, Parra T, Aubourg C, Hébert R, Philippot P. Strain localization during crustalscale boudinage t ‑ o form extensional metamorphic domes in the Aegean Sea. Special papers—Geological Society of America; 2004. p. 185–210. Joly A, Martelet G, Chen Y, Faure M. A multidisciplinary study of a syntectonic pluton close to a major lithospheric‑scale fault—Relationships between the Montmarault granitic massif and the Sillon Houiller Fault in the Variscan French Massif Central: 2. Gravity, aeromagnetic investigations, and 3‑D geologic modeling. J Geophys Res. 2008;113:1. Király Á, Capitanio F‑A, Funiciello F, Faccenna C. Subduction induced mantle flow: length‑scales and orientation of the toroidal cell. Earth Planet Sci Lett. 2017;479:284–97. Kushnir A‑R, Heap M ‑ J, Baud P. Assessing the role of fractures on the permeability of the Permo‑ Triassic sandstones at the Soultz‑sous‑Forêts (France) geothermal site. Geothermics. 2018;74:181–9. Kranz R‑L, Frankel A‑D, Engelder T, Scholz C‑H. The permeability of whole and jointed Barre granite. Int J Rock Mech Min‑ ing Sci Geomech. 1979;16(4):225–34. Lahsen A, Muñoz N, Parada M‑A. Geothermal development in Chile. In: Proceedings world geothermal congress, Bali, Indonesia. 2010. Launay G. Hydrodynamique des systèmes minéralisés péri‑ granitiques: étude du gisement à W‑Sn‑(Cu) de Panasqueira (Portugal). Dissertation, University of Orléans. 2018. Lebocey J. Géologie et gîtologie du district du Pontgibaud (Puy‑ de‑Dôme). In: Bayle LD, De Ascençao Guedes R, Gol D, editors. Le Règne Minéral‑hors série XIX ‑2013. p 35–42. Ledésert B. Fracturation et paléocirculations hydrothermales. Application au granite de Soultz‑sous‑Forêts, Dissertation, Université de Poitiers. 1993. López D‑L, Smith L. Fluid flow in fault zones: analysis of the interplay of convective circulation and topographically driven groundwater flow. Water Resour Res. 1995;31(6):1489–503. Lowell R‑P. Circulation in fractures, hot springs, and convective heat transport on mid‑ ocean ridge crests. Geophys J Int. 1975;40(3):351–65. Lucazeau F. Flux de chaleur, production de chaleur et évolution géodynamique récente du Massif Central Français. Dis‑ sertation, Université des sciences et techniques du Languedoc, Montpellier. 1981. Lucazeau F, Vasseur G, Bayer R. Interpretation of heat flow data in the French Massif Central. Tectonophysics. 1984;103(1):99–119. Lucazeau F, Vasseur G. Heat flow density data from France and surrounding margins. Tectonophysics. 1989;164(2):251–8. Magri F, Möller S, Inbaretal N. 2D and 3D coexisting modes of thermal convection in fractured hydrothermal systems— implications for transboundary flow in the Lower Yarmouk Gorge. Mar Pet Geol. 2016;78:750–8. Malkovsky V‑I, Magri F. Thermal convection of temperature ‑ dependent viscous fluids within three ‑ dimensional faulted geothermal systems: estimation from linear and numerical analyses. Water Resour Res. 2016;52:2855–67. Malkovsky V‑I, Pek A‑A. Timescales for reaching steady‑state fluid flow in systems perturbed by formation of highly permeable faults. Geofluids. 2004;4:253–8. Manning C‑E, Ingebritsen S‑E. Permeability of the continental crust: implications of geothermal data and metamorphic systems. Rev Geophys. 1999;37(1):127–50. Mareschal J‑ C, Jaupart C. Energy budget of the Earth. In: Gupta HK (ed) Encyclopedia of solid earth geophysics. Dordrecht: Springer; 2011. p. 285–91. Martel C, Champallier R, Prouteau G, Pichavant M, Arbaret L, Balcone‑Boissard H, Scaillet B. Trachyte phase rela‑ tions and implication for magma storage conditions in the Chaîne des Puys (French Massif Central). J Petrol. 2013;54(6):1071–107. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 29 of 29 Mazurek M. Geological and hydraulic properties of water‑ conducting features in crystalline rocks. Hydrogeology of crystalline rocks. Dordrecht: Springer; 2000. p. 3–26. Mckenna J‑R, Blackwell D ‑D. Numerical modeling of transient basin and range extensional geothermal systems. Geother ‑ mics. 2004;33(4):457–76. Meredith P‑ G, Main I‑ G, Clint O‑ C, Li L. On the threshold of flow in a tight natural rock. Geophys Res Lett. 2012;39:4. Moeck I‑ J. Catalog of geothermal play types based on geologic controls. Renew Sustain Energy Rev. 2014;37:867–82. Mordensky S‑P, Villeneuve M ‑ C, Kennedy B‑M, Heap M ‑ J, Gravley D‑M, Farquharson J‑I, Reuschlé T. Physical and mechani‑ cal property relationships of a shallow intrusion and volcanic host rock, Pinnacle Ridge, Mt. Ruapehu, New Zealand. J Volcanol Geothermal Res. 2018;359:1–20. Micarelli L, Benedicto A, Wibberley CAJ. Structural evolution and permeability of normal fault zones in highly porous carbonate rocks. J Struct Geol. 2006;28(7):1214–27. Murphy H‑D. Convective instabilities in vertical fractures and faults. J Geophys Res. 1979;84:6121–30. Négroni J‑M. Le district de Pontgibaud cadre Géologique Evolution Structurale et Métallogénique. Dissertation, Univer ‑ sité de Clermont‑Ferrand. 1981. Neuman S‑P. Generalized scaling of permeabilities: validation and effect of support scale. Geophys Res Lett. 1994;21(5):349–52. Oda M. An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses. Water Resour Res. 1986;22(13):1845–56. Person M, Hofstra A, Sweetkind D, Stone W, Cohen D, Gable C‑ W, Banerjee A. Analytical and numerical models of hydro‑ thermal fluid flow at fault intersections. Geofluids. 2012;12(4):312–26. Pribnow D, Schellschmidt R. Thermal tracking of upper crustal fluid flow in the Rhine Graben. Geophys Res Lett. 2000;27(13):1957–60. Rabinowicz M, Boulègue J, Genthon P. Two and three‑ dimensional modeling of hydrothermal convection in the sedi‑ mented middle valley segment, Juan de Fuca ridge. J Geophys Res. 1998;103(B10):24045–65. Ranjram M, Gleeson T, Luijendijk E. Is the permeability of crystalline rock in the shallow crust related to depth, lithology or tectonic setting? Geofluids. 2015;15(1–2):106–19. Roche V, Sternai P, Guillou‑Frottier L, Menant A, Jolivet L, Bouchot V, Gerya T. Emplacement of metamorphic core com‑ plexes and associated geothermal systems controlled by slab dynamics. Earth Planet Sci Lett. 2018a;49:322–33. Rowland J‑ V, Sibson R‑H. Structural controls on hydrothermal flow in a segmented rift system, Taupo Volcanic Zone, New Zealand. Geofluids. 2004;4(4):259–83. Roche V, Bouchot V, Beccaletto L, Jolivet L, Guillou‑Frottier L, Tuduri J, Tokay B. Structural, lithological, and geodynamic controls on geothermal activity in the Menderes geothermal Province ( Western Anatolia, Turkey). Int J Earth Sci. 2018b;104:1–28. Sausse J. Caractérisation et modélisation des écoulements fluides en mi‑ lieu fissuré. Relation avec les altérations hydro ‑ thermales et quantification des paléo ‑ contraintes. Dissertation, Université Henri Poincaré ‑Nancy I. 1998. Sardini P, Ledésert B, Touchard G. Quantification of microscopic porous networks by image analysis and measurements of permeability in the Soultz‑sous‑Forêts granite (Alsace, France). In: Jamtveit B, Yardley BWD, editors. Fluid flow and transport in rocks Mechanisms and effects. New York: Chapman & Hall; 1997. p. 171–89. Smith L, Chapman DS. On the ther ‑ mal effects of groundwater flow: 1. regional scale systems. J Geophys Res. 1983;88(1):593–608. Siebenaller L, Boiron M‑ C, Vanderhaeghe O, Hibsch C, Jessel M, André‑Mayer A‑S, Photiades A. Fluid record of rock exhumation across the brittle‑ ductile transition during formation of a Metamorphic Core Complex (Naxos Island, Cyclades, Greece). J Metamorph Geol. 2013;31:313–38. Sobolev S, Zeyen H, Granet M, Achauer U, Bauer C, Werling F, Altherr R, Fuschs F. Upper mantle temperatures and lithosphere‑asthenosphere system beneath the French Massif Central constrained by seismic, gravity, petrologic and thermal observations. Tectonophysics. 1997;275(1–3):143–64. Sonney R, Vuataz F‑D. Numerical modelling of Alpine deep flow systems: a management and prediction tool for an exploited geothermal reservoir (Lavey‑les‑Bains, Switzerland). Hydrogeol J. 2009;17(3):601–16. Süer S, Güleç N, Mutlu H, Hilton D‑R, Çifter C, Sayin M. Geochemical monitoring of geothermal waters (2002–2004) along the North Anatolian Fault Zone, Turkey: spatial and temporal variations and relationship to seismic activity. Pure Appl Geophys. 2008;165(1):17–43. Taillefer A, Guillou‑Frottier L, Soliva R, Magri F, Lopez S, Courrioux G, Le Goff E. Topographic and faults control of hydro ‑ thermal circulation along dormant faults in an orogen. Geochem Geophys Geosyst. 2018;19:4972–95. Taillefer A, Soliva R, Guillou‑Frottier L, Le Goff E, Martin G, Seranne M. Fault ‑related controls on upward hydrothermal flow; an integrated geological study of the Têt fault, eastern Pyrénées (France). Geofluids. 2017;2017:19. Vasseur G, Gable R, Feuga B, Bienfait G. Groundwater flow and heat flow in an area of mineral springs. Geothermics. 1991;20(3):99–117. Verney N. Importance du contexte géologique dans la signature chimique et isotopique (Carbone 13) des sources hydrothermales. Applications aux émergences du Puy de Dôme. Université Blaise Pascal, MsC thesis. 2005. Violay M, Heap M‑ J, Acosta M, Madonna C. Porosity evolution at the brittle‑ ductile transition in the continental crust: implication for deep hydro‑ geothermal circulation. Sci Rep. 2017;7(1):7705. Vigneresse J‑L, Tikoff B, Améglio L. Modification of the regional stress field by magma intrusion and formation of tabular granitic plutons. Tectonophysics. 1999;302(3–4):203–24. Walter B. Réservoir de socle en contexte extensif: genèse, géométries et circulations de fluids: exemples du rift intraconti‑ nental du lac Albert (Ouganda) et de la marge proximale d’Ifni (Maroc). Dissertation, Université de Lorraine. 2016. Zeyen H, Novak O, Landes M, Prodehl C, Driad L, Hirn A. Refraction‑seismic investigations of the northern Massif Central (France). Tectonophysics. 1997;275(1–3):99–117. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geothermal Energy Springer Journals

On the geothermal potential of crustal fault zones: a case study from the Pontgibaud area (French Massif Central, France)

Loading next page...
 
/lp/springer-journals/on-the-geothermal-potential-of-crustal-fault-zones-a-case-study-from-FRPkZ35S4L
Publisher
Springer Journals
Copyright
Copyright © 2019 by The Author(s)
Subject
Earth Sciences; Geotechnical Engineering & Applied Earth Sciences; Renewable and Green Energy; Geoecology/Natural Processes
eISSN
2195-9706
DOI
10.1186/s40517-019-0150-7
Publisher site
See Article on Publisher Site

Abstract

hugo.duwiquet@gmail.com TLS‑ Geothermics, 92 The present study aims to understand the potential of a new and novel type of geo‑ Chemin de Gabardie, thermal play system for high temperature and electricity production: crustal fault zones 31200 Toulouse, France Full list of author information (CFZ). According to geological and geophysical data, the Pontgibaud fault zone (French is available at the end of the Massif Central) is suspected to host an active hydrothermal system at a depth of a few article kilometers. The deep geometry of the fault zone and the permeability distribution are the main unknown parameters that are required to assess the geothermal potential of the Pontgibaud site. Structural and thin‑section observations, laboratory permeability and connected porosity measurements and X‑ray micro ‑tomography observations suggest that the hydrothermal system behaves like a double matrix‑fracture perme ‑ ability reservoir. Numerical modeling in which we varied the fault dip and the ratio between the fault zone permeability and host rock, R, was performed. Results indicate that three main convective regimes can be identified (weak convection, single cellular ‑ type convection and bicellular convection). For a sufficiently high fault zone permeabil‑ −15 2 ity (> 1 × 10 m ), buoyancy‑ driven flow creates a positive thermal anomaly of several tens of °C at a depth of 2–5 km. For a vertical fault zone, the thermal anomaly is larger for higher R values. Numerical models, then applied to the geologically constrained Pontgibaud fault zone, show that a temperature of 150 °C at a depth of 2500 m can be −14 2 obtained for a fault zone permeability of 1.6 × 10 m . Based on a multi‑ disciplinary approach, this work establishes a potential predictive tool for future high‑temperature geothermal operations within basement rocks hosting large‑scale fault systems. Keywords: Crustal fault zones, Permeability, High‑temperature geothermal system, Numerical modeling, Pontgibaud, French‑Massif Central Introduction Exploring new geothermal targets requires the understanding of the factors affecting fluid circulation and heat transfer (Rowland and Sibson 2004; Fairley 2009). Understand - ing the distribution of permeability in the crust remains an essential component for the general comprehension of a geothermal model in the crustal domain, and therefore for the success of a geothermal prospect (Moeck 2014). Nowadays in France, most geother- mal exploration licences are located in sedimentary basin areas or within graben struc- tures. The Paris and Aquitaine basins, characterized by simple sedimentary structures, allow geothermal exploitation from low (30 °C) to medium (110 °C) temperature systems. © The Author(s) 2019. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 2 of 29 Geothermal exploitation is also rife in the northeast of France in the Upper Rhine Graben. The enhanced geothermal systems at Soultz-sous-Forêts and Rittershoffen (Alsace, France) reach temperatures up to about 160  °C (e.g., Pribnow and Schellschmidt 2000; Baujard et al. 2017). However, no exploitation is currently being developed in the basement domain that host crustal-scale faults that are likely to modify fluid circulation. Numerous studies highlight that crustal-scale fault zones located in an extensional regime represent efficient conduits for meteoric fluids to flow down to mid-crustal depths (Barton et al. 1995; Faulds et al. 2010; Mckenna and Blackwell 2004; Siebenaller et al. 2013; Haines et al. 2016; Roche et al. 2018a, b). At these depths (8–15 km), corresponding to, or close to, the brittle–duc- tile transition, temperatures may exceed 400 °C (Violay et al. 2017). Fault zones appear as areas of high permeability in which fluids can circulate and focus temperature anomalies at the surface (Faulds et al. 2010; Taillefer et al. 2018a, b). However, the generation of a geo- thermal resource in a fault zone likely depends on its geometry with respect to the present- day stress field, on its permeability structure and on the local thermal regime. The present study aims at understanding the potential of a new and novel type of geothermal system for electricity production: the crustal fault zone (CFZ). Crustal fault zones can be defined as crustal-scale heterogeneities, which localize the deformation and thus modify the mechanical properties of the crust, down to the brit- tle–ductile transition (BDT), whose depth depends on the local thermo-mechanical regime. In some cases, CFZ may cross BDT and extend below it (e.g., Famin et al. 2004; Jolivet et  al. 2004). CFZ should not be considered as single major fault but rather as a large deformation zone (several hundreds to thousands of meter wide). These crustal fault zones can be found worldwide and their geothermal potential can be highlighted by the occurrence of geothermal events, e.g., North Anatolian Fault Zone, Turkey (Süer et al. 2008), Liquine-Ofqui fault zone, Chile (Lahsen et al. 2010), Badenweiler-Lenzkirch Sutur, Germany (Brockamp et al. 2015). However, their geothermal potential has never been studied and even less quantified. Among the first studies on fluid circulation within a fault, Murphy (1979) proposed a theoretical approach using the Rayleigh number. He highlighted vertical thermal insta- bilities linked to convection cells. Convective heat transfer establishes positive and nega- tive thermal anomalies within a fault (Lowell 1975; Murphy 1979; Rabinowicz et al. 1998; Faulds et al. 2010). While many parameters seem to condition a convective regime, the permeability ratio between the fault and its host rock emerges as a key parameter (For- ster and Smith 1989; López and Smith 1995). While fault zones often host hydrothermal systems (Faulds et al. 2010), the permeability of these networks depends on a multitude of parameters. The relationships between (i) fault geometry and stress field (Barton et al. 1995); (ii) fractures and their connectivity, and (iii) alteration, dissolution and precipita- tion, are all processes that can cause faults to behave as conduits for (Caine et al. 1996; Mazurek 2000; Gudmundsson 2000; Grasby and Hutcheon 2001; Brogi and Fulignati, 2012; Belgrano et  al. 2016) or barriers to (Bense et  al. 2008; Bertrand 2016) fluid flow. Consequently, it is worth using a multidisciplinary approach to strongly constrain fluid flow models in hydrothermal systems (Oda 1986; Caine et  al. 1996; Bense and Person, 2006; Micarelli et  al. 2006; Faulkner et  al. 2010; Bense et  al. 2013). As an example, the recent study by Edel et al. (2018) predicts, using a combination of magnetic, gravity, seis- mic and geological data, a temperature of 150 °C at a depth of 2500 m in the northern Duwiquet et al. Geotherm Energy (2019) 7:33 Page 3 of 29 Vosges area (France). With an improved numerical modeling approach, Guillou-Frottier et al. (2013) reproduced and predicted temperatures at geothermal sites of Soultz-sous- Forêts and Rittershoffen. Since 2014, TLS-Geothermics, a French company involved in geothermal exploration has been keen to demonstrate the viability of CFZ as a geothermal exploration play for economic power generation. The demonstration of the potential of CFZ is underway at a site in the French Massif Central (FMC) using an exploration licence named “La Sioule”. Several geological observations and geophysical data, which are detailed below, represent good indicators for high-temperature geothermal exploitation (tempera- tures > 150 °C for shallower than 3500 m). The Pontgibaud fault zone, located in the La Sioule licence (Fig.  1) corresponds to a well-suited case study since numerous data are available (structural geology, lithology, topography, springs geochemistry, thermal properties, magnetotellurics, gravimetry, ambient noise tomography, and seismicity) (Bellanger 2017). Thanks to this complete dataset, and with new additional petrophysical data, 2D numerical modeling can help to better understand the links between fluid flow and permeability. The possible presence Fig. 1 Simplified geological map showing the “La Sioule” licence area (black contour) and a brief overview of the major stratigraphic units in the study area. French Massif Central (FMC). Coordinate system: WGS84 Pseudo‑Mercator EPSG:3857 (After the 1/50,000th geological map published by BRGM) Duwiquet et al. Geotherm Energy (2019) 7:33 Page 4 of 29 of convection cells is here examined to understand the establishment of positive or nega- tive thermal anomalies in the fault system. Permeability will be considered as a spatially variable parameter to analyze its effect on the hydrothermal regime of the fault system. Geological setting Geological formations The La Sioule licence area (black contour in Fig.  1) is located to the east of the Limagne graben in the Massif Central, in a Paleozoic basement domain and partially covered by the northern part of the Cenozoic volcanic formations (in blue in Fig.  1). It is bounded to the west by the Sillon Houiller fault zone, to the south by the volcanic massif of the Monts Dore, and to the north and east by the Carboniferous basin of Manzat and the granitic massif of Saint-Gervais, respectively. Variscan metamorphic formations (ortho and para-migmatites, orthogneiss, paragneiss and mica schists) are mainly present in the area (in green in Fig. 1). Carboniferous plutonic rocks are characterized by monzo- granite and granodiorite (red in Fig.  1) and Viseo-Serpukhovian volcano–sedimentary formations (brown in Fig.  1). Upper Pennsylvanian and Permian detritic sedimen- tary formation mark the late Variscan basins emplaced on the major crustal fault zone (mapped by the BRGM). Structures and favorable features for fluid circulation in the “La Sioule” licence The Sillon Houiller fault zone This fault zone, generally oriented N 20°, is located to the west of the area (Fig.  2) and is continuous at a scale of 400  km with a sinistral motion (Faure et  al. 1993, Joly et  al. 2008). Evidence in various volcanic deposits highlights the high permeability of this structure since the Permian. At both sides of the Sillon Houiller, P-wave velocity anoma- lies indicate that this structure affect the entire lithosphere (Granet et  al. 1995a). Such a large-scale structure, if permeable enough, would promote the ascent of hot hydro- thermal fluids from the base of crust to the surface. According to Sobolev et al. (1997), temperature varies between the different sides of the fault at 30  km depth, from 700 −2 to 925  °C. In addition, surface heat flow equals 113 mW  m near the Sillon Houiller (Fig. 2) (Lucazeau et al. 1984; Lucazeau and Vasseur 1989). This abnormally high value is thought to be related to a relatively high mantle heat flow (Lucazeau et al. 1984), which could be related to a supposed mantle plume (Granet et al. 1995b) or past subduction- related hot upwelling (Király et al. 2017; Roche et al. 2018ab). The Aigueperse–Saint‑Sauves fault zone This NE–SW fault system is more than 300 km long (Faure et al. 1993) and crosses the southeast corner of the licence (Fig.  2). It comprises a fault network, several hundred meters thick, and shifts the Sillon Houiller fault zone over a 31-km-long dextral-slip motion. The Puy de la Nugère and Puy Chopine monogenic volcanoes of the quater - nary Chaîne des Puys (CdP) volcanic field (located 10 km to the east of La Sioule licence area) are emplaced on the Aigueperse–Saint-Sauves fault zone. The presence of Ceyssat, Duwiquet et al. Geotherm Energy (2019) 7:33 Page 5 of 29 Fig. 2 Thermal, hydrothermal and geological features in the studied area. Location of Pranal and Ceyssat springs (brown dots), drilling of Peyrouses 1 (black dot), heat flow and geothermal gradient (red crosses, Lucazeau and Vasseur 1989; Vasseur et al. 1991; International Heat Flow Commission). Gelles granite (Négroni 1981). Heat production (Lucazeau 1981) and resistivity profile (black dashed line, Ars, 2018) Châtelguyon, Gimaux, Saint Myon or Montpensier springs characterizes this structure as a potentially open and efficient pathway for hydrothermal fluids. The Pontgibaud fault zone The N–S fault zone of Pontgibaud is parallel to the CdP. The CdP is a fissural volcanic field that hosts about eighty basaltic to trachytic volcanoes that erupted between 100 and 7 ka. This system (Fig.  2), more than 40 km long, is pinched between the Sillon Houiller and the Aigueperse–St-Sauve fault zones. In the studied area, it bounds the Viseo-Ser- pukhovian Gelles granite (Fernandez 1969; Négroni 1981). Some N–S faults of the Pon- tgibaud fault zone are filled by microgranite. Moreover, this N–S network is affected by As–Au, Pb–Ag and Ba–F mineralization, reflecting the ability of the system to transmit hydrothermal fluids during the Variscan. Some dating on mineralized bodies obtained Duwiquet et al. Geotherm Energy (2019) 7:33 Page 6 of 29 on phengite and illite suggest Permian to Jurassic ages (Bril et  al. 1991), which reflects late hydrothermal fluid circulation related to the opening of the Alpine Tethys and the Atlantic Ocean. This mineralization is attributed to a medium temperature phase (200– 400 °C) sub-contemporary to the end of the crystallization of the Gelles granite. If evi- dence for hydrothermal fluid circulation exists during the Carboniferous, but also from the Permian to Jurassic, current fluid circulation along this fault zone seems also very likely. Indeed, springs can be found along the present-day fault at several locations (see details below). Geothermal setting of the Pontgibaud area Springs and magma chambers Present-day fluid flow occurs within the Pontgibaud fault zone (Fig.  2). In addition, at the time of mining activities, the galleries were regularly flooded (Lebocey 2013). Due to cessation of mining activities in the 20th century, the mining wells were abandoned. Some wells, such as the one at Sainte-Barbe in Barbecot (4  km north of Pontgibaud), are marked by ion-rich and extremely CO -charged flows. In addition, the Pranal spring (Fig.  2) comes from a deep reservoir with estimated temperatures between 150 and 175 °C (estimated using Na/K geothermometer; Verney 2005). At a deeper level, thermo-barometric experiments on the CdP trachytes estimated the storage pressure and temperature conditions of magmas (Martel et  al. 2013). For temperatures ranging from 700 to 825  °C, the results of Martel et  al. (2013) indicate pressures between 3 and 3.5 kbar (i.e., between 11 and 13  km deep, respectively). The thermal relaxation time calculated for these reservoirs indicate that cooling started about 15,000  years ago (Martel et  al. 2013). This means that they are still hot (500– 600 °C), if not partially melted. Surface heat flow and heat production Part of the Pontgibaud fault zone was intersected by core drilling undertaken by the BRGM (French Geological Survey) in the 1970s. The purpose of these drillings was to characterize the mining potential of these structures. As part of this study, one of these boreholes (Peyrouses 1, location in Fig.  2) remains accessible and has been used in this study. In the past, temperature measurements have been carried out between 55 and 120 m. From these measurements, a geothermal gradient of 41 °C/km has been inferred −2 and a heat flux of 110  mW  m was estimated (IHFC database; http://www.datap ages. c om/g i s -ma p- publi s hing - pr o g r am/g i s -op en-file s /g lob a l-f rame wor k/g lob a l-he a t -flow - datab ase). This relatively high heat flow value can be explained by an elevated mantle heat flow and/or a large crustal heat production. The analysis of 660 rock samples taken from granitic units in the FMC (Lucazeau 1981) revealed an average heat production of −3 4.06 μW m . Considering an average granite thickness of 7 km (Vigneresse et al. 1999), the heat production of this granitic upper crust can only explain 26% of the heat flow measured in the Pontgibaud fault zone. The analysis of 112 samples taken from meta - −3 morphic units revealed an average heat production of 2.27 μW m . If this average value is representative for the remaining 22 km of crust (considering the crustal thickness in Duwiquet et al. Geotherm Energy (2019) 7:33 Page 7 of 29 −2 this area is 29 km, Zeyen et al. 1997), the mantle heat flow would be close to 32 mW m , twice as large as for a thermally stable Archean crust. Permeability In a crustal context, permeability and connected porosity values are generally considered −18 very low to medium, with porosity ranging from 0 to 15% and permeability from 10 −13 2 to 10  m (Brace 1984; Evans et  al. 1997; Sonney and Vuataz 2009; Ingebritsen and Appold 2012; Ranjram et  al. 2015; Walter 2016; Achtziger-Zupancic et  al. 2016, 2017). However, the orders of magnitude for permeability values and derived parameters vary greatly from one site to another and need to be addressed specifically to better charac - terize a specific site or case study (Manning and Ingebritsen 1999). Therefore, new per - meability measurements on representative samples of the various lithologies from the Pontgibaud fault zone have been performed. An additional way to constrain the fault zone system is to use geophysical data, as detailed below. Magnetotelluric data Magnetotelluric (MT) campaigns were carried out in the Pontgibaud fault zone area by TLS-Geothermics and IMAGIR between 2015 and 2017 (Ars et al. 2019). From their 3D model of electrical resistivity produced in 2016 by IMAGIR, the profile corresponding to our geological cross-section was extracted (Fig.  3). The interest of this passive elec - tromagnetic method is its sensitivity to geometries and its ability to image conductive objects in three dimensions. Results show that the obtained conductive bodies in the upper crust could be related to the Pontgibaud fault zone. Mining work on the Pontgibaud fault zone has highlighted structures with a dip of 70° E (Bouladon et al. 1961), but more recent data favors dips of ~ 60°, as detailed later. This Pontgibaud fault zone Aigueperse-Saint-Sauves fault zone Fig. 3 2D resistivity profile. The dashed lines represent the hypothetical tracing of the deep structures, constrained by surface signatures (fault orientations and dips). (Ars et al. 2019) Duwiquet et al. Geotherm Energy (2019) 7:33 Page 8 of 29 steep dip is consistent with the structures that are visible on the MT image (Fig. 3). The obtained resistivity model makes it possible to identify highly conductive bodies, one of which is located between a depth of 10 and 20 km (Fig. 3). Petrological observations on trachytes suggest that the conductivity anomaly located at a depth between 10 and 20 km can reflect the presence of silicate melt (Martel et al. 2013). At a depth of 8 and 3  km, two other anomalies located in the brittle crust could be linked to the presence of magmas, clays (smectites) and/or hydrothermal fluids. What - ever their nature, these two anomalies appear to be related to the Pontgibaud fault zone. This fault zone, which could be interpreted as a listric fault zone, seems to join the Aigueperse–Saint-Sauves fault system at a depth of 10 to 15  km (Ars et  al. 2019). Isotopic analyses of CO from the Ceyssat and Pranal resurgences at δ C (Fig.  2) indi- cate a mantle origin (Verney 2005), which is in agreement with isotopical results from the whole French Massif Central (Brauer et al. 2017). Indeed, the presence of CO at the surface can be explained by the connectivity of the Pontgibaud and Aigueperse–Saint- Sauves structures to magmatic activity. It also highlights that these two fault zones serve as pathways for hydrothermal fluid flow. Methods Field observations This approach aims to characterize the intensity of the brittle deformation as well as the general geometry of the Pontgibaud fault zone. Structural data (fault orientations and dips) collected since 2014 by TLS-Geothermics have been compiled to assess the inten- sity of deformation. These observations are key to improve the geometry and scaling of the simplified numerical model. The microstructural analysis of the samples collected within this hydrothermal zone (“Peyrouses 1” borehole). Figure  2 provides an under- standing of how fluids flow within this deformed zone. Laboratory observations and measurements In order to qualitatively and quantitatively characterize the permeability and porosity parameters of the hydrothermal zone, 2D thin-section observations, 3D X-ray micro- tomography observations and laboratory measurements of permeability and connected porosity were performed. The thin-section observations were performed on the same samples that were analyzed by X-ray micro-tomography. X-ray micro-tomography is a non-destructive imaging technique that provides a three-dimensional visualization of the analyzed samples. X-rays emitted by a tungsten filament pass through the sample. Its power depends on the thickness and nature of the radiographed element. A detector, located behind the sample, measures the attenua- tion of the X-ray beam. A series of images is then generated in different greyscales. A three-dimensional reconstruction of the object is performed using this series of images. However, the density contrast that may exist between dense mineral phases (pyrite) and low-density mineral phases (feldspar) can bias the information. Beam–Hardening–Cor- rection (BHC) limits this loss of information. This imaging technique makes it possible to obtain voxel resolutions between 25 and 4  μm. Once rebuilt, the three-dimensional volume analysis and processing can be performed using the Volume Graphics VG Studio Max software (https ://www.volum egrap hics.com/en/produ cts/vgstu dio-max.html). Duwiquet et al. Geotherm Energy (2019) 7:33 Page 9 of 29 The connected porosities and permeabilities of selected samples were measured at the Ecole et Observatoire des Sciences de la Terre (EOST) in Strasbourg (France). Measure- ments were performed on the seven samples for which the porosity was visualized using X-ray micro-tomography (including fractured samples, intact samples, altered samples, and altered and fractured samples). First, the oven-dried (in a vacuum at 40 °C for at least 48 h) 20 mm diameter and 40 mm long core samples were placed in an AccuPyc II helium pyc- nometer to calculate their skeletal (connected volume). The connected porosity (the inter - connected porosity within the sample) was then calculated using the bulk volume of the samples, measured using digital calipers. Second, the gas permeability was measured using a benchtop nitrogen permeameter (Farquharson et al. 2016; Heap and Kennedy 2016; Heap 2019). All measurements were made at room temperature and under a confining pressure of 1 MPa. Permeability was measured using either the “steady-state” or the “pulse-decay” method, depending on the permeability of the sample [methods are outlined in detail in Heap et al. (2017) and Kushnir et al. (2018)]. For the steady-state measurements, volumetric flow rates were measured (using a gas flow meter) for different values of differential pressure (upstream minus the downstream pressure; measured using a pressure transducer). For the pulse-decay measurements, the pressure decay of a known volume of gas (in the upstream circuit) was monitored as a func- tion of time. Estimate of Rayleigh number The Rayleigh number expresses the ratio between the driving forces (e.g., gravity) that pro - mote fluid circulation and the resistant forces (e.g., viscous resistance) that prevent circula - tion. It is expressed in the form: KLα�Tg Ra = , (1) Dν 2 −1 where K (m ) is the medium permeability, L (m) is the thickness of the system, α (°C ) is the thermal expansion coefficient, ∆T (°C) is the temperature difference between the −2 upper and lower limits, g is the acceleration due to gravity (m s ) and ν is the kinematic 2 −1 2 −1 fluid viscosity (m  s ). D is the thermodispersion tensor (m  s ) (e.g., Magri et al. 2016) as defined by: Φf + (1 − Φ)s D = , (2) ρf Cpf −1 −1 where Φ corresponds to the porosity, Cpf (J  kg  K ) to the fluid thermal capacity, λf −1 −1 and λs, respectively, to the fluid and solid thermal conductivity (W m  K ). The critical Rayleigh number is the value above which heat is evacuated through convec - tive circulation. For 3D vertical faults with a permeable top, and in the case of a constant fluid viscosity (not temperature-dependent), Malkovsky and Pek (2004) provided an analyt - ical relationship to calculate the critical Rayleigh number (Rcf ). This relationship depends only on the fault aspect ratio Δ: 0.8584 1.165 6.428 1.165 Rcf = + (27.1) , (3) � Duwiquet et al. Geotherm Energy (2019) 7:33 Page 10 of 29 where , where d is the thickness of the fault zone (m) and H its depth (m). Other 2H numerical and analytical simulations have shown that, for a temperature-dependent fluid viscosity, the critical permeability above which convection begins is lower than the case where the fluid viscosity is considered constant (Malkovsky and Magri 2016). According to the calculations of Malkovsky and Magri (2016), the critical fault perme- ability Kfc from which convection is likely to be initiated in the faults is calculated using Rcf /4. Numerical approach In order to quantify the effect of fault zone dip and the impact of variable permeability on fluid circulation in the fault zone, a series of numerical simulations using Comsol Multiphysics software was performed. Heat equations and Darcy’s law were coupled in a similar way to the Guillou-Frottier et al. (2013). The numerical approach was organized in two successive stages. The first stage was dedicated to the study of a theoretical fault zone and focused on the influence of fault dip and permeability on fluid flow in order to understand the processes governing deep fluid circulation. The second stage was dedi - cated to the Pontgibaud hydrothermal system, in order to verify the legitimacy of the digital crustal-scale model. Finally, the obtained results, and in particular the geothermal gradient and surface heat flow, were compared with field data, including the resistivity data from Ars et  al. (2019) (Fig.  3) which were not considered to build our numerical model. Boundary and initial conditions The average temperature at the surface was fixed at 10  °C. According to regional data (see “Geothermal setting of the Pontgibaud area” section), a constant heat flux of 35 −2 mW  m was imposed at the base of the model. The lateral borders of the model were thermally insulated. The initial conditions of the model were defined from the stationary state obtained for the pure thermal conduction regime. At t = 0, the permeabilities were assigned to the different compartments, and the fluid flow was initiated. Considering the thermal relaxation time of the magmatic chambers (see above, Martel et  al. 2013) calculations were performed up to a time of t + 15,000 years. For Darcy’s law, a no flow condition was imposed at the bottom and on both sides of the model, a surface pressure equivalent to atmospheric pressure (P = 10 Pa) was imposed at the surface and a positive pressure evolution with depth was adopted as the initial condition. Meshing and geometry Two geometrical configurations were tested. The first configuration was that of the theo - retical fault zone (Fig.  4a), and the second was the crustal-scale model applied to the Pontgibaud case study (Fig.  4b), which is constrained by observations and field meas - urements. For the smallest (i.e., for the models of the theoretical fault zone) fault zone models, the mesh consisted of 61,211 triangles with a size not larger than 5  m in the high-permeability zones and not smaller than 150 m in the low-permeability zones. The mesh was also refined at the intersections of the fault planes. For the largest model (i.e., Duwiquet et al. Geotherm Energy (2019) 7:33 Page 11 of 29 Fig. 4 Geometries studied by numerical modeling. a Theoretical fault zone geometry. b Simplified crustal‑scale geometry containing essential structural elements derived from field and geophysical data, and MT interpretation the crustal models), the mesh consisted of 52,798 triangles with a size not larger than 25 m in the high-permeability zones and not smaller than 300 m in the low-permeability zones. Variable physical properties The dynamic viscosity of the fluid is function of temperature T (°C), such that: −5 µ (T ) = 2.414 · 10 · exp (4) T + 133 according to Smith and Chapman (1983) and Rabinowicz et  al. (1998). The density of water is also a function of temperature: ρ (T ) = 1002.4 − 0.1905 · T − 0.0025 · T (5) based on a polynomial adjustment of experimental data up to 350 °C, as used by Guil- lou-Frottier et  al. (2013). For temperatures above 350  °C, the fluid density was fixed at −3 600  kg  m . This hypothesis can be considered valid since such temperatures within Duwiquet et al. Geotherm Energy (2019) 7:33 Page 12 of 29 −16 2 the fault are reached at depths where permeability is lower than 10 m , preventing buoyancy forces to be effective. The permeabilities K of the basement (m ) and K of the b f faults (m ) are the main adjustment variables. For the basement, the decrease in perme- ability with depth can be physically characterized by the following equation (Garibaldi et al. 2010): z − 800 K (z) = K exp , (6) b b where K (z) is the permeability of the depth-dependent basement, K is the permeabil- b b ity at the surface, and z is the depth below sea level. The length δ (m) characterizes the intensity of the decrease in permeability with depth, taken here at 2500 m. As the Gelles granite and the pre-Tertiary bedrock are intensely fractured, Eq. (6) was applied to these lithologies, in accordance with the models of Garibaldi et al. (2010) and Guillou-Frottier et  al. (2013). According to laboratory measurements (see below), a permeability value −16 2 of 10 m was assigned to K . In addition, to account for the uncertainty associated with upscaling the experimental values (Neuman 1994; Heap and Kennedy 2016), a large range of values has been tested in the numerical models. Within the fault zone, we considered that the permeability variation with depth varied according to Eq.  (7). However, field observations indicate that the area has a very high fracture density. In addition, the lateral variation in permeability will be very heterogeneous. Therefore, we consider that the use of Eq.  (6) is not suffi- cient. By keeping the decrease in permeability with depth, a lateral variation can be included, by using the following spatial variation: z − 800 z + x ( ) ( ) K (x, z) = K × exp × 101 + 100 × sin 2 × π × , f f (7) where K (x, z) is the space-dependent permeability of the fault and K × 201 being the f f maximum permeability value at the surface. The last term on the right side of Eq.  (7) makes the permeability alternate along a sinusoid applied to the vertical and horizontal axes of the numerical model. The term λ corresponds to the wavelength of the sinusoid. To reproduce rather fine alternations of high and low permeability, as suggested by field observations (Fig.  6), a low value of λ was chosen. Finally, the choice of constants 100 and 101 makes it possible to vary the permeability over two orders of magnitude. Numerical modeling allows us to test the effect of the unknown parameters. In order to see the effects of geometry and permeability, the physical properties on either side of the fault zone were identical and fixed (Fig. 5). During the various sim- ulations, the dip of the fault zone was tested at angles ranging from 10° to 90°. The assigned permeability to the host and fault zone varied according to Eqs. (6) and (7). −18 −16 2 The K value varied from 1 × 10 to 2 × 10 m , thus giving a maximum perme- −14 2 ability value ( K ) of 4 × 10 m , hence greater than the critical permeability esti- max −16 2 mated for the Pontgibaud fault zone. For the basement, the K value was 10 m −20 2 and decreased down to 1 × 10 m at depth. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 13 of 29 Fig. 5 Permeability imposed in the numerical models of the parametric study. The alpha angle corresponds to the dip of the fault zone (here α = 90°). The permeability of the basement and fault will decrease with depth according to Eq. 6 (Garibaldi et al. 2010). Within the fault zone, permeability will also vary laterally (Eq. 7) to match field observations Results Structure characterization With 126 new measurements of fault dips and orientations, three families of faults have been identified. The first family comprises 67° ESE dipping normal N 10° trending vein beam faults (Fig.  6a). With a surface thickness of 3 km, this normal fault zone has low- ered the topography by 300  m eastwards. In addition, over the study area, the Sioule River surprisingly flows north–south over more than 10 km, and partly borders the vein network, before becoming meandering again. This north–south direction is also consist - ent with the main direction of the Limagne fault, which is located a little further east. Due to field data and geophysical arguments (see “Magnetotelluric data ” section), the general geometry of the Pontgibaud fault zone could be regarded as a listric type. The second and third fault families, both steep, intersect the Pontgibaud fault zone. The second family is formed of mainly dextral N 30° trending faults (Fig.  6b), while the third family consists of sinistral faults that trend N 115° (Fig. 6c). To summarize all field observations and data from the literature, a compiled geologi - cal cross-section is shown in Fig.  7. Its location corresponds to the geophysical profile (Fig. 3) and it crosses the villages of Prondines, Ceyssat, and ends at the Limagne basin. This cross-section illustrates all the geological and geophysical features detailed in the previous sections. Thus, the proposed listric geometry of Pontgibaud fault zone is asso - ciated with two other vertical fault families. The presence of springs within these struc - tures is an indicator of present-day fluid circulation. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 14 of 29 Fig. 6 Observations of the different lithologies and fault families present in the study area. a′, c′ Metamorphic series. b′ Granite. Schmidt’s projection on the lower hemisphere (a–c). Isodensity obtained with Stereonet ( TLS‑ Geothermics, modified) Fig. 7 Synthetic and simplified geological cross‑section of the Pontgibaud area Laboratory observations and analyses Qualitative approach: thin‑section observations Thin-section observations show fracture-rich and highly altered facies (Fig.  8). These fractures are mainly filled with quartz, feldspar, phyllosilicates, ferriferous, pyrite and oxide minerals. If at the heart of these fractures the circulation of fluids appears diffi - cult, the walls seem to be favorable for the circulation of fluids (Fig.  8c). Fluid circulation within the fractures promoted the crystallization of secondary mineral phases. Turma- line (Fig.  8d), is observed in large quantities in all thin-sections. Fluid circulation also caused significant mineral dissolution (Fig.  8b). This dissolution can be accompanied by a loss of volume, thus promoting the formation of pores within the matrix (Launay 2018). These initial observations show that fractures related to tectonic history served as pathways for (paleo)fluids and produced haloes of alteration. These haloes were subject to significant mineral dissolution resulting in the formation of micro-porosity. Fractures and micro-porosity represent the available space for modern-day fluid circulation. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 15 of 29 Fig. 8 Microscopic observations of thin‑section from the hydrothermalized zone. a Anatexite rich in veins. b Fractured anatexite. c Mineralized breach. d Porous breach. (Po) Porosity. (Qz) Quartz. ( Void) Void. (Fr) Fracture. ( Tur) Turmaline. (Ox) Oxide. (Sr) Sericites. (Kfs) Potassic feldspar However, the ability of fluids to flow through the hydrothermal system will depend on the connectivity of these fractures and micro-porosity. Three-dimensional imaging is required to qualitatively estimate the connectivity of these microstructural elements. Qualitative approach: X‑ray micro‑tomography The results of the reconstruction and image processing are given in Fig.  9. Figure  9a shows the full sample with a voxel resolution of 25  μm. It is characterized by a matrix (grey in Fig.  9a) and by fractures (yellow and red in Fig.  9a). Two phases are present within these fractures. The red phase characterizes the mineralization of the fracture. Thin-section observations on this sample showed that this fracture is filled with lead sulphides. The yellow phase represents the planar porosity (i.e., fractures). The pro - cessing of this image (Fig.  9a) made it possible to isolate these fractures (Fig.  9b) and qualitatively characterize the connectivity of the planar porosity. The three-dimensional visualization (Fig.  9b) shows a high density of fractures. These fractures are distributed heterogeneously in space. This three-dimensional arrangement promotes connectivity of the pores located at the fracture walls (yellow in Fig.  9b). Because a visualization of the matrix porosity (Fig.  9b) was not possible at this resolution, due to the characteristic pore size, a smaller volume of this sample was analyzed with a resolution of 4 μm in an attempt to detect and segment the potentially connected pores. The result of the analysis and processing of the sample is given in Fig. 9c. Coupling thin-section (2D) and micro- tomography (3D) observations highlight altered planar (i.e., fractures) and matrix (i.e., pores) porosities partially filled with secondary minerals. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 16 of 29 Fig. 9 Three‑ dimensional observations of the planar and matrix porosity of one sample taken from the Pontgibaud fault zone. The red color of figures a, b corresponds to mineralization (lead sulphides). The yellow color of figures a, b corresponds to the planar porosity located on the fracture walls. The yellow color of figure c could correspond to the matrix porosity. The resolution of figures a–c is, respectively, 25, 25 and 4 μm Because this small-scale porous network could play a crucial role in determining the fluid circulation within the Pontgibaud hydrothermal zone, we performed direct meas - urements of connected porosity and permeability in the laboratory. Quantitative approach: porosity and permeability measurements −18 −13 2 The measured permeability values varied between 2 × 10 and 7 × 10 m while the connected porosity values varied between 6 and 22% (Fig. 10). In general, permeability increases as connected porosity increases (Fig. 10). The high - est porosity and permeability were measured for the fractured and altered sample, and the lowest porosity and permeability were measured for the intact and unaltered sample (Fig. 10). To summarize, the three above-described results suggest that permeability is increased by the presence of system-lengthscale fractures and alteration (Fig.  10). The relatively high permeability of the measured samples (the permeability of granite can be as low −22 2 as 10 m ; Meredith et  al. 2012), suggest that fluid flow within the matrix could also be significant, alongside large-scale fractures and fracture networks. Hence, after hav - ing determined the first unknown parameter in “Structure characterization ” section Duwiquet et al. Geotherm Energy (2019) 7:33 Page 17 of 29 Fig. 10 Connected porosity and permeability measurements on sample taken from the Peyrouse 1 borehole (identical to those used for qualitative observations). The non‑fractured and unaltered samples have the lowest permeability and connected porosity value. Conversely, the altered and fractured sample has the highest value of permeability and connected porosity. These purely qualitative degrees of alteration and fracturation are based on macroscopic, microscopic observations and three‑ dimensional visualization by X‑ray micro ‑tomography analysis (fault geometry), we identify here the second one (permeability variation). To better understand how fault geometry and permeability variation influence hydrothermal fluid flow and thus control the establishment of thermal anomalies, a numerical approach is required. Numerical approach Rayleigh number analysis Using the Pontgibaud fault zone parameters (H = 15,000 m and h = 3000 m) and Eq. (3), the critical Rayleigh number Rcf /4 equals 21.3. The measured permeability varies from −18 −13 2 10 to 10 m (Fig.  10), implying a variation of the Rayleigh number (Eq.  1) from 0.018 to 1800. We can thus deduce the critical fault permeability above which thermal convection can occur: Rcf · D · ν −15 2 (8) K = = 1.2 × 10 m , fc g · α · �T · L −4 −1 ◦ −6 2 −1 −6 2 −1 where α = 2 × 10 K , T = 450 C , D = 0.76 × 10 m s , and ν = 10 m s . u Th s, considering the selected parameters, it seems reasonable to consider a convective thermal regime within the Pontgibaud fault zone, at least in the case when the perme- −15 2 ability is greater than 1.2 × 10 m . This threshold value must, however, be considered with caution since the fluid properties considered here correspond to those for room temperature conditions. Numerical simulations of a theoretical fault zone More than 200 simulations were performed, with the numerical set-up of Fig.  5. The results are summarized in Fig. 11. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 18 of 29 Fig. 11 Graphical representation of all numerical models of temperature anomalies on the theoretical fault zone performed (n = 216). The values of positive thermal anomalies are shown in red. Their depths are represented in brown The images in Fig.  11 represent temperature anomalies (∆T), corresponding to abnormally hot and abnormally cold temperatures compared to an undisturbed (con- ductive) thermal regime. The amplitude of the temperature anomalies is different for each case, but numerical values are indicated for some of them. In order to account for both the permeability of the fault and of the basement, the R ratio is defined by: max R = . (9) This ratio compares the maximum permeability of the fault and of the basement −16 2 (the permeability K being set at 10 m ). Other K values have been tested, but b b 0 0 higher values are not realistic and lower values do not change the results detailed below. As it can be seen in Fig. 11, for R = 200, as the dip of the fault zone increases from 10° to 90°, the temperature values increase from 12 to 85.4 °C (red numbers in Fig. 11). Thus, for a fixed R ratio, the highest thermal anomaly value will be for a vertical fault zone. The depth of the thermal anomaly also varies as a function of fault dip. By increasing the fault dip from 10° to 90°, the depth of the thermal anomaly (brown numbers in Fig. 11) decreases from 3.6 to 0.7 km. Thus, for a fixed R value, vertical structures will result in the largest thermal anomalies at the shallowest depths. By observing the variations in the temperature anomalies, we can see that when R increases for a fixed angle of 60°, the values of the thermal anomalies are 68.6  °C, 82.5 and 77.4  °C. Indeed, an increase in the R ratio does not show a significant increase in the temperature anomaly. However, the depth of these same anomalies decreases from 2.7 to 0.5 km. So for a fixed dip, when R ratio is high, the thermal anomaly will be at a maximum at a shallower depth. Therefore, a fault zone with a large dip and a high R ratio will result in the largest thermal anomalies at the shallowest depths. While this repre- sentation characterizes the distribution and intensity of positive thermal anomalies, it does not characterize the processes at their origin. Based on the results of these same numerical simulations, Fig.  12 graphically summarizes the morphology of the numeri- cally obtained isotherms as a function of fault dip angle, value of positive thermal anom- aly and the R ratio. Three main zones can be described: Duwiquet et al. Geotherm Energy (2019) 7:33 Page 19 of 29 Fig. 12 Diagram of convective regimes. Grey lines represent isotherm morphology. The dashed black line separate three area. In blue, the area of unicellular weak type convection zone. In orange, the area of unicellular medium‑type convection zone. In red, the area of bicellular strong‑type convection zone. The dashed grey contour on the right of the diagram illustrates the possible range for the Pontgibaud fault zone (i) For R values between 1.9 and 250 (blue area in Fig.  12) the morphology of iso- therms is slightly deformed. This low deformation occurs at dip values below 50°. The value of the resultant thermal anomaly will be lower than 65 °C. These results are consistent with Fig.  11 since the lower the dip of the structure and the lower the R ratio, the lower the thermal anomaly produced. The black arrows show that the circulation of fluids is controlled by a single long-wavelength convection cell. This domain will be called the “unicellular weak-type convection zone”. (ii) For R values between 40 and 400 with a dip value between 10 and 90° (orange area in Fig.  12), the isotherms show moderate deformation. The range of values of the thermal anomalies produced is between 10 and 120 °C. In the same way, the arrows show that fluid circulation is controlled by moderate-wavelength convection cells. This domain will be called the “unicellular medium-type convection zone”. (iii) Finally, for R values between 120 and 400 and dips greater than 50° (red area in Fig. 12) the isotherms are significantly deformed. The value of the thermal anomaly is higher than 67  °C. These results also show consistency with Fig.  11 since the higher the dip and R ratio, the greater the thermal anomaly at shallow depth pro- duced. The black arrows show that the circulation of fluids is controlled by two small-wavelength convection cells. This domain will be called the “bicellular strong-type convection zone”. The results of this parametric study should be applicable to hydrothermal systems in a crustal context. To sum up, for unicellular convection, the values of the thermal anomalies produced are between 10 and 120 °C. For bicellular convection, the values of the thermal anomalies produced are between 70 and 120  °C. In the next section, we apply these results to the case study of the hydrothermal system at Pontgibaud. According to Fig.  12, the Pontgibaud fault zone should fall somewhere within these two domains. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 20 of 29 Numerical simulations of the Pontgibaud fault zone The geometry of the numerical model for Pontgibaud (Fig.  13) is directly based on the geological model (Fig.  5) and is therefore constrained by a detailed set of geological data and observations. The physical properties, integrated into the numerical models, are summarized in Appendix A. In accordance with the thermal relaxation time of the magmatic chambers calculated on the basis of a thermo-barometric study (Martel et al. 2013), the results of our numerical simulations are shown up to time t + 15,000 years. Figure  13 presents the results of the numerical simulations with the temperature field (colors) and isotherms in white. The stable temperature pattern shown in Fig.  13a −15 2 is obtained for K value up to 1.4 × 10 m . For these K values, the isotherms f f max max are not deformed, but are disturbed around the cooling magma chambers. At time Fig. 13 Maps of temperatures (a), permeabilities (b), isotherms (c), and thermal anomalies (d) of the Pontgibaud region (at time t + 15,000 years). In accordance with Fig. 11, color of boxes contours correspond to one convective regime: blues boxes correspond to areas of unicellular weak‑type convection zone, oranges boxes correspond to unicellular medium‑type convection zone, and reds boxes correspond to bicellular strong‑type convection zone Duwiquet et al. Geotherm Energy (2019) 7:33 Page 21 of 29 t + 15,000 years, the cooling of the magma chambers does not disturb any isotherm in the Pontgibaud fault zone. Figure  13b represents the permeability field centered on the red box of Fig. 13a. The deformed isotherms shown in Fig.  13c, d are located at the inter- section between the Pontgibaud listric geometry and the vertical faults (see black box, −15 −14 2 Fig.  13a). For K values between 1.6 × 10 and 1.6 × 10 m , the deformation of max the isotherms can be observed (Fig.  13c). In the terms previously used, the morphol- ogy of the isotherms changes from areas of weak convection of unicellular type (blue −15 2 box, K = 1.6 × 10 m ) to areas of strong convection of bicellular type (red box, max −14 2 K = 1.6 × 10  m ). This representation therefore shows that, when permeability is max increased within the Pontgibaud fault zone, isotherms show an increasing disturbance. Figure  13d shows the amplitudes and locations of the obtained thermal anomalies. These anomalies are represented as a function of R ratio, the value of the maximum tem - perature anomaly (∆T ) and its depth. As it can be seen in Fig. 13d, positive and nega- max tive thermal anomalies are present in the Pontgibaud system. Positive thermal anomalies are located at the top of the fault zone, while negative thermal anomalies are located at the base of the fault zone. The depth of the positive thermal anomaly decreases for an increasing R ratio (see dashed line in Fig.  13d). This latter observation is consistent with the parametric study (Figs.  11, 12). For R ratios between 42 and 160, the value of the thermal anomaly does not show significant variation. Overall, these anomalies are between 57 and 66 °C. On the other hand, when R = 16, the value of the thermal anom- aly is only 24 °C. This lower value can be explained by a low deformation of the isotherms. At R = 42 and above, the isotherms show an increasing deformation up to R = 160. Considering isotherms with a relatively large deformation, the value of the positive thermal anom- aly no longer varies significantly. Thus, the results of numerical simulations carried out within the Pontgibaud hydrothermal system show that the deformation of the isotherms −15 2 is significant from a K value of 1.6 × 10 m (in accordance with the previous Ray- max leigh number analysis) and a R ratio of 16. Combining field data, laboratory measurements and hydrothermal models Finding the right R and K combination for Pontgibaud fault zone The comparison of the results obtained by modeling with the data measured in the field makes it possible to assess the consistency of the obtained results. Within the Pontgibaud −2 fault zone, the surface heat flux is 110 mW m and the geothermal gradient is 41 °C/km (Vasseur et al. 1991, International Heat Flow Commission). This comparison will make it possible to estimate the permeability value for which the field measurements and simu - lation results are consistent. Among the numerous tested models, it turns out that only some combinations precisely reproduced the above-described thermal features. The −14 2 result for the simulation with a K = 1.6 × 10 m , show a surface heat flux of 115 max −2 −1 −14 2 mW m and a geothermal gradient of 39 °C km . Thus, for K = 1.6 × 10 m , the max results of the numerical models reproduce the heat flux and geothermal gradient values measured near the surface. The numerical models show that at a time t + 10,000  years −2 the heat flux and the geothermal gradient are 105 mW m and 33 °C/km, respectively. −2 At t + 20,000 years we obtained 127 mW m and 41 °C/km, respectively. 0 Duwiquet et al. Geotherm Energy (2019) 7:33 Page 22 of 29 Combining geophysical data and numerical results Figure  14a, b and c compare the temperature anomaly map of the Pontgibaud region with the measured resistivity profile extracted at the exact coordinates of the geological cross-section. The map (Fig.  14a) shows positive thermal anomalies within the Pontgibaud fault zone, at the magma chambers and at the base of the sedimentary infill of the Limagne basin (right part of the model). A negative thermal anomaly is also present at the base of the Pontgibaud fault zone. The MT image shows various anomalies of low resistivity at 3, 8 and 15 km depth (Fig.  14b). By superimposing the two images (Fig.  14a, b), we can see that the temperature anomalies located at 3, 8 and 15 km depth are superimposed on the low-resistivity anomalies located at 3, 8 and 15 km depth. At 15 km depth, the thermobarometric study by Martel et al. (2013) shows the pres- ence of a magma chamber. At 3 and 8  km depth, Fig.  14c shows that anomalies of low resistivity are located in areas where fluids could circulate. Thus, the numerical simula - −14 2 tion imposing a permeability K = 1.6 × 10 m at the Pontgibaud fault zone shows max a striking consistency with the resistivity profile, even if the link between resistivity and Fig. 14 a Thermal anomalies maps of Pontgibaud area. The black arrows represent the direction of fluid flow. b Resistivity model previously described (Fig. 3). c Comparison of the numerical simulation result affecting −14 −2 a K permeability of 1.6 × 10 m at the Pontgibaud fault zone, with the resistivity profile. d Fluid flow max velocity within the Pontgibaud fault zone Duwiquet et al. Geotherm Energy (2019) 7:33 Page 23 of 29 Fig. 15 Numerical simulation results for a maximum permeability imposed on the Pontgibaud fault zone of −14 −2 1.6 × 10 m temperature may not be so obvious. Nevertheless, isotherms deformation is correlated with high fluid flow velocity (Fig. 14d). Discussion Multi‑disciplinary approach This study aims to understand how crustal-scale fault zones can generate a viable geo - thermal resource. The overall understanding of an outcropping fault system inte - grates many parameters that can be studied using a multi-disciplinary approach. Thus, structural observations and measurements, 2D and 3D thin-section and X-ray micro- tomography observations, as well as connected porosity and permeability laboratory measurements, have highlighted two essential variables: the fault dip angle and perme- ability of a fault zone. Numerical modeling has made it possible to test the effect of these two variables on the intensity and depth of positive temperature anomalies, but also to characterize the deformation of the resultant isotherms. The comparison of the results of the numerical simulations at crustal scale with the field data shows a permeability −14 2 of the Pontgibaud fault zone K = 1.6 × 10 m . For this value, the ratio R = 160 max (Figs.  11, 12). The value of the thermal anomaly varies between 68 and 70  °C and has depths between 2.7 and 1.7  km (Figs.  11, 12). Our results show a thermal anomaly of 66 °C within the Pontgibaud fault zone at a depth of 1.8 km (Fig.  15). It is worth men- −14 2 tioning that the permeability value of 1.6 × 10 m is in accordance with previous Duwiquet et al. Geotherm Energy (2019) 7:33 Page 24 of 29 studies where temperature measurements and surface heat flow values were numerically reproduced with similar permeability values (e.g.,  Person et  al. 2012; Guillou-Frottier et al. 2013; Roche et al. 2018a, b). Origin of porosity Thin-sections observations (Fig.  8) revealed open fractures and fractures completely sealed by mineral precipitation. A fault can thus behave as a pathway for, or as a barrier to, fluid circulation (Bense et  al. 2013). Quartz, feldspar and alteration products such as sericites are well known to form around hydrothermal sites and near faults (Bruhn et  al. 1994; Sausse 1998; Caine et  al. 2010). Filling fractures with cement or alteration products can reduce their permeability (e.g., Griffiths et al. 2016; Mordensky et al. 2018). When fractures transmit fluids, their immediate surroundings can be subject to signifi - cant mineral dissolution (potassic feldspar, biotite, plagioclases), forming matrix poros- ity around these fractures. If the minerals likely to be altered are initially connected in space, then their alterations will produce connected pores that can promote the circula- tion of fluids (e.g., Farquharson et al. 2018). Direct connected porosity measurements on our samples with no apparent fractures (at the lengthscale of the 40  mm-long sample) showed values between 6 and 22%. Indeed, the permeability of the altered samples was higher than that of the unaltered samples (Fig.  10). The porosity values measured here are consistent with the study by Sardini et al. (1997) who quantified the connectivity of the minerals likely to be altered. Mercury porosity between 6 and 20% was measured on samples of altered, non-fractured granite. Considering the capacity of fluids to circulate through fractures but also within the matrix itself, it is thus possible to consider Pon- tgibaud hydrothermal system as a reservoir with double matrix-fracture permeability. This kind of reservoir has already been observed in the altered granite of the Soultz- Sous-Forêts geothermal reservoir (Genter 1989; Ledésert 1993). Controlling fluid flow factors within the Pontgibaud hydrothermal system Although we did not consider the present-day stress field, the value of thermal anom - alies and their depths are the result of convective heat transfer processes. López and Smith (1995) showed that this thermal regime is directly related to the permeability con- trast between the fault and its host rocks. In addition to the results of López and Smith (1995), in which the R ratio was studied, we have explored here the effect of the dip angle of the fault zone on fluid circulation. Putative structures sufficiently permeable to allow fluid circulation are commonly described for hydrothermal systems (Forster and Evans 1991; Bruhn et al. 1994; Belgrano et al. 2016) and are used in many numerical models (Forster and Smith 1989; Mckenna and Blackwell 2004; Garibaldi et al. 2010; Guillou-Frottier et al. 2013; Magri et al. 2016; Taillefer et al. 2018a, b). The results of the parametric study presented herein, and their comparison with the Pontgibaud natural system, show that when the fault zone is verti- cal, the value of thermal anomaly will be high and at shallow depth. This depth is also controlled by the ratio between the fault and its host permeability: when the R ratio is high, thermal anomaly will be at a shallower depth. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 25 of 29 Permeability measurements were carried out on samples from the Pontgibaud fault −18 −13 2 zone (Peyrouses 1 drilling). The results show values between 2 × 10 and 7 × 10 m . The permeabilities of the rocks measured in this study are relatively high compared to laboratory measurements on intact granite samples (e.g., Brace et al. 1968; Kranz et al. 1979; Meredith et al. 2012). The maximum permeability value for the basement assigned −16 2 to the ratio R has been set at 10 m . This value is compatible with the permeabil - ity values extracted from the laws of permeability decrease proposed by Ingebritsen and Appold (2010, 2012). Conclusions The aim of this study was to better understand the organization and controlling factors of a hydrothermal system associated with crustal faults, with an application to the Pon- tgibaud hydrothermal system. Based on results from field observations, microscopic observations and laboratory analyses, the Pontgibaud fault zone has the characteristics of a dual matrix-fracture permeability reservoir. The results also highlighted parameters that are essential for understanding a hydrothermal system: the dip angle of the fault and its permeability. Numerical 2D simulations, combining fluid flow and the heat equation, have tested the effect of these parameters on the value and distribution of positive ther - mal anomalies within the fault zone. Based on more than 200 numerical simulations, the results of this parametric study showed that (i) a vertical structure provides the largest thermal anomaly at the shallowest depth; (ii) the depth of positive thermal anomaly is shallow for high R ratio values. Although the model is quite simple, it well captures the convective behavior. These generic results were then compared to a natural case study, the hydrothermal system at Pontgibaud. Numerical simulations performed at a crustal scale, and integrat- ing observations and field measurements, showed similar results to those of the afore - mentioned parametric study. However, a 3D study should bring new insights into the Pontgibaud hydrothermal system and its associated thermal anomalies. The results of fluid circulation simulations within the Pontgibaud fault zone showed −14 2 that a maximum fault permeability value of 1.6 × 10 m allows reproduction of the surface data. Considering the economic potential, a temperature of 150 °C at a depth of 2500  m represents an interesting target for high-temperature geothermal exploitation. Finally, a new geothermal play, crustal fault zones was identified and intensively explored and analyzed. A methodology has been developed to constrain this specific hydrother - mal system, and to target economically interesting areas. Furthermore drilling planned by operators should make it possible to verify and validate all or parts of this approach. Abbreviations BRGM: Bureau de Recherches Géologiques et Minières; CFZ: Crustal fault zone; EOST: Ecole et Observatoire des Sciences de la Terre; FMC: French Massif Central; ISTO: Institut des Sciences de la Terre d’Orléans; IMAGIR: https ://www.imagi r.eu/; TLS: TLS‑ Geothermics (https ://www.tls‑geoth ermic s.fr). Acknowledgements The financial support was provided by TLS‑ Geothermics and supported by ISTO and BRGM. Funding for open access was generously offered by the Helmholtz Centre for Environmental Research‑UFZ, Helmholtz Centre Potsdam‑ GFZ German Research Centre for Geosciences, the Karlsruhe Institute of Technology (KIT ), and LABEX Grant ANR‑11‑LABX ‑0050_G‑ EAU‑ THERMIE‑PROFONDE (this research therefore benefited from state funding managed by the Agence National de la Recherche (ANR) as part of the “Investissements d’avenir” program)”. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 26 of 29 Authors’ contributions This work corresponds to M.Sc. thesis of HD. LA analyzed, interpreted and retrieved the field data. LGF contributed to numerical modeling and participated in writing the manuscript. MJH helped HD to perform the permeability measure‑ ments. MB provided a very large amount of data and participated in the various field campaigns. All authors read and approved the final manuscript. Funding The financial support was provided by TLS‑ Geothermics. Availability of data and materials All data are available in Figs. 1, 2, 4, 8 and Appendix A. Competing interests The authors declare that they have no competing interests. Author details 1 2 3 Université d’Orléans, ISTO, UMR 7327, 45071 Orléans, France. BRGM, 45060 Orléans, France. ISTO, UMR7327, Université d’Orléans, CNRS, BRGM, 45071 Orléans, France. TLS‑ Geothermics, 92 Chemin de Gabardie, 31200 Toulouse, France. Géophysique Expérimentale, Institut de Physique de Globe de Strasbourg, UMR 7516, CNRS, Université de Strasbourg/ EOST, 5 Rue René Descartes, 67084 Strasbourg Cedex, France. Appendix Appendix A Rock properties Symbols Magmatic Calcareous‑ Metamorphic Gelles Carboniferous Migmatitic– Mantle Unit chamber clay unit granite plutonic rocks gneissic complex complex Porosity Ф 1.0 10 5 5 5 3 2 % Perme‑ k Variable Variable Variable Variable Variable Variable Vari‑ m ability able Thermal λ 1.3 2.2 2.5 3.25 2.8 2.5 3 W/m K con‑ duc‑ tivity Heat Cp 800 880 800 800 800 800 1000 J/(kg K) capac‑ ity −3 Heat A 0.5 1.7 2.1 3.5 3.5 2 0.5 µW m pro‑ duc‑ tion Bulk ρ 3300 2650 2700 2600 2600 3000 3300 kg/m den‑ sity Fluid properties Symbols Values Unit Bulk density ρ Variables kg/m Thermal conductivity λ 0.6 W/m K Heat capacity Cp 4200 J/(kg K) Dynamic viscosity µ Variables Pa s Properties of these parameters have been recovered in the following works: Mareschal and Jaupart (2011), Lucazeau et al. (1984), Hasterok et al. (2017), Taillefer (2017), Guil- lou-Frottier et  al. (2013), Launay (2018), McKenna and Blackwell (2004), Smith and Chapman (1998), Rabinowicz et al. (1998). Duwiquet et al. Geotherm Energy (2019) 7:33 Page 27 of 29 Received: 21 May 2019 Accepted: 24 October 2019 References Achtziger‑Zupančič P, Loew S, Hiller A, Mariethoz G. 3D fluid flow in fault zones of crystalline basement rocks (Poehla‑ Tellerhaeuser Ore Field, Ore Mountains, Germany). Geofluids. 2016;16(4):688–710. Achtziger‑Zupančič P, Loew S, Mariéthoz G. A new global database to improve predictions of permeability distribution in crystalline rocks at site scale. J Geophys Res. 2017;122(5):3513–39. Ars JM, Tarits P, Hautot S, Bellanger M, Coutant O, Maia M. Joint inversion of gravity and surface wave data constrained by magnetotelluric: application to deep geothermal exploration of crustal fault zone in felsic basement. Geothermics. 2019;80:56–68. Barton C‑A, Zoback M ‑D, Moos D. Fluid flow along potentially active faults in crystalline rock. Geology. 1995;23(8):683–6. Baujard C, Genter A, Dalmais E, Maurer V, Hehn R, Rosillette R, Schmittbuhl J. Hydrothermal characterization of wells GRT‑1 and GRT ‑2 in Rittershoffen, France: implications on the understanding of natural flow systems in the rhine graben. Geothermics. 2017;65:255–68. Belgrano T‑M, Herwegh M, Berger A. Inherited structural controls on fault geometry, architecture and hydrothermal activ‑ ity: an example from grimsel pass, Switzerland. Swiss J Geosci. 2016;109(3):345–64. Bellanger M. High temperature geothermal resources of crustal fault zones: a dedicated approach. In: 79th EAGE confer‑ ence and exhibition, June 2017. https ://doi.org/10.3997/2214‑4609.20170 1771. Bense V, Person M, Chaudhary K, You E, Cremer Y, Simon N‑S. Thermal anomalies indicate preferential flow along faults in unconsolidated sedimentary aquifers. Geophys Res Lett. 2008;35(24):1025. Bense V‑F, Person M ‑A. Faults as conduit ‑barrier systems to fluid flow in siliciclastic sedimentary aquifers. Water Resour Res. 2006;42:5. Bense V, Gleeson T, Loveless S, Bour O, Scibek J. Fault zone hydrogeology. Earth Sci Rev. 2013;127:171–92. Bertrand L. Etude des réservoirs géothermiques développés dans le socle et à l’interface avec les formations sédimen‑ taires. Dissertation, Université de Lorraine. 2016. Bouladon JP, Picot J‑ J, Sainfeld P. Le faisceau filonien de Pontgibaud. Bull BRGM. 1961;1:1–41. Brace A‑ W. Permeability of crystalline rocks: new in situ measurements. J Geophys Res. 1984;89(B6):4327–30. Brace W, Walsh J‑B, Frangos W ‑ T. Permeability of granite under high pressure. J Geophys Res. 1968;73(6):2225–36. Bril H, Bonhomme M‑ G, Marcoux E, Baubron J‑ C. Ages K/Ar des minéralisations de Brioude‑Massiac ( W ‑Au‑As‑Sb; Pb ‑Zn), Pontgibaud (Pb‑Ag; Sn), et Labessette (As‑Pb ‑Sb ‑Au): place de ces districts dans l’évolution géotectonique du Massif central français. Miner Deposita. 1991;26(3):189–98. Brockamp O, Schlegel A, Wemmer K. Complex hydrothermal alteration and illite K‑Ar ages in Upper Visean molasse sedi‑ ments and magmatic rocks of the Variscan Badenweiler‑Lenzkirch suture zone, Black Forest, Germany. Int J Earth Sci. 2015;104(3):683–702. Brogi E, Fulignati P. Tectonic control on hydrothermal circulation and fluid evolution in the pietratonda‑poggio peloso (southern Tuscany, Italy) carbonate‑hosted sb ‑mineralization. Ore Geol Rev. 2012;44:158–71. Bruhn R‑L, Parry W ‑ T, Yonkee W‑A, Thompson T. Fracturing and hydrothermal alteration in normal fault zones. Pure Appl Geophys. 1994;142(3–4):609–44. Caine J‑S, Evans J‑P, Forster C‑B. Fault zone architecture and permeability structure. Geology. 1996;24(11):1025–8. Caine J‑S, Bruhn R‑L, Forster C‑B. Internal structure, fault rocks, and inferences regarding deformation, fluid flow, and mineralization in the seismogenic stillwater normal fault, dixie valley, nevada. J Struct Geol. 2010;32(11):1576–89. Edel J‑B, Maurer V, Dalmais E, Genter A, Richard A, Letourneau O, Hehn R. Structure and nature of the Palaeozoic base ‑ ment based on magnetic, gravimetric and seismic investigations in the central Upper Rhinegraben. Geotherm Energy. 2018;6(1):13. Evans J‑P, Forster C‑B, Goddard J‑ V. Permeability of fault‑related rocks, and implications for hydraulic structure of fault zones. J Struct Geol. 1997;19(11):1393–404. Fairley J‑P. Modeling fluid flow in a heterogeneous, fault ‑ controlled hydrothermal system. Geofluids. 2009;9(2):153–66. Famin V, Philippot P, Jolivet L, Agard P. Evolution of hydrothermal regime along a crustal shear zone, Tinos Island, Greece. Tectonics. 2004. https ://doi.org/10.1029/2003T C0015 09. Farquharson J‑I, Wild B, Kushnir A‑R, Heap M ‑ J, Baud P, Kennedy B. Acid‑induced dissolution of andesite: evolution of permeability and strength. J Geophys Res. 2018. https ://doi.org/10.1029/2018J B0161 30. Farquharson J‑I, Heap M ‑ J, Lavallée Y, Varley N‑R, Baud P. Evidence for the development of permeability anisotropy in lava domes and volcanic conduits. J Volcanol Geoth Res. 2016;323:163–85. Faulds J, Coolbaugh M, Bouchot V, Moek I, Oguz K. Characterizing structural controls of geothermal reservoirs in the Great Basin, USA, and western turkey: Developing successful exploration strategies in extended terranes. In: Pro‑ ceedings world geothermal congress Bali, Indonesia, April 25–30. 2010. Faulkner D, Jackson C, Lunn R, Schlische R, Shipton Z, Wibberley C, Withjack M. A review of recent developments con‑ cerning the structure, mechanics and fluid flow properties of fault zones. J Struct Geol. 2010;32(11):1557–75. Faure M, Grolier J, Pons J. Extensional ductile tectonics of the Sioule metamorphic series ( Variscan French Massif Central). Geol Rundsch. 1993;82(3):461–74. Fernandez A. La série cristallophyllienne et les granites de la région de Pontgibaud (Puy de Dôme) Massif Central Fran‑ çais. Dissertation, Université de Clermont Ferrand. 1969. Forster C, Smith L. The influence of groundwater flow on thermal regimes in mountainous terrain: a model study. J Geophys Res. 1989;94(B7):9439–51. Forster C‑B, Evans J‑P. Hydrogeology of thrust faults and crystalline thrust sheets: results of combined field and modeling studies. Geophys Res Lett. 1991;18(5):979–82. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 28 of 29 Garibaldi C, Guillou‑Frottier L, Lardeaux J‑M, Bonté D, Lopez S, Bouchot V, Ledru P. Thermal anomalies and geologi‑ cal structures in the Provence basin: implications for hydrothermal circulations at depth. Bull Soc Géol France. 2010;181(4):363–76. Genter A. Géothermie Roches Chaudes Sèches: le granite de Soultz‑sous‑Forêts. (Bas‑Rhin, France). Dissertation, Univer ‑ sité d’Orléans. 1989. Granet M, Stoll G, Dorel J, Achauer U, Poupinet G, Fuchs K. The Massif Central (France). New constraints on the geody‑ namical evolution from teleseismic tomography. Geophys J Int. 1995a;121:33–48. Granet M, Wilson M, Achauer U. Imaging a mantle plume beneath the French Massif Central. Earth Planet Sci Lett. 1995b;136:281–96. Grasby S‑E, Hutcheon I. Controls on the distribution of thermal springs in the southern cordillera. Can J Earth Sci. 2001;38(3):427–40. Griffiths L, Heap M ‑ J, Wang F, Daval D, Gilg H‑A, Baud P, Genter A. Geothermal implications for fracture ‑filling hydrother ‑ mal precipitation. Geothermics. 2016;64:235–45. Gudmundsson A. Active fault zones and groundwater flow. Geophys Res Lett. 2000;27(18):2993–6. Guillou‑Frottier L, Carre C, Bourgine B, Bouchot V, Genter A. Structure of hydrothermal convection in the upper Rhine graben as inferred from corrected temperature data and basin‑scale numerical models. J Volcanol Geoth Res. 2013;256:29–49. Hasterok D, Webb J. On the radiogenic heat production of igneous rocks. Geosci Front.2017;8(5):919–40. Heap M‑ J. The influence of sample geometry on the permeability of a porous sandstone. Geosci Instrum Methods Data Syst. 2019;8(1):55–61. Heap M‑ J, Kennedy B‑M. Exploring the scale ‑ dependent permeability of fractured andesite. Earth Planet Sci Lett. 2016;447:139–50. Heap M‑ J, Violay M, Wadsworth F‑B, Vasseur G. From rock to magma and back again: the evolution of temperature and deformation mechanism in conduit margin zones. Earth Planet Sci Lett. 2017;463:92–100. Haines S, Lynch E, Mulch A, Valley J‑ W, van der Pluijm B. Meteoric fluid infiltration in crustal‑scale normal fault systems as indicated by δ18O and δ2H geochemistry and 40Ar/39Ar dating of neoformed clays in brittle fault rocks. Litho‑ sphere. 2016;8(6):587–600. Ingebritsen S‑E, Manning C‑E. Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism. Geofluids. 2010;10(1–2):193–205. Ingebritsen S‑E, Appold M ‑S. The physical hydrogeology of ore deposits. Econ Geol. 2012;107(4):559–84. Jolivet L, Famin V, Mehl C, Parra T, Aubourg C, Hébert R, Philippot P. Strain localization during crustalscale boudinage t ‑ o form extensional metamorphic domes in the Aegean Sea. Special papers—Geological Society of America; 2004. p. 185–210. Joly A, Martelet G, Chen Y, Faure M. A multidisciplinary study of a syntectonic pluton close to a major lithospheric‑scale fault—Relationships between the Montmarault granitic massif and the Sillon Houiller Fault in the Variscan French Massif Central: 2. Gravity, aeromagnetic investigations, and 3‑D geologic modeling. J Geophys Res. 2008;113:1. Király Á, Capitanio F‑A, Funiciello F, Faccenna C. Subduction induced mantle flow: length‑scales and orientation of the toroidal cell. Earth Planet Sci Lett. 2017;479:284–97. Kushnir A‑R, Heap M ‑ J, Baud P. Assessing the role of fractures on the permeability of the Permo‑ Triassic sandstones at the Soultz‑sous‑Forêts (France) geothermal site. Geothermics. 2018;74:181–9. Kranz R‑L, Frankel A‑D, Engelder T, Scholz C‑H. The permeability of whole and jointed Barre granite. Int J Rock Mech Min‑ ing Sci Geomech. 1979;16(4):225–34. Lahsen A, Muñoz N, Parada M‑A. Geothermal development in Chile. In: Proceedings world geothermal congress, Bali, Indonesia. 2010. Launay G. Hydrodynamique des systèmes minéralisés péri‑ granitiques: étude du gisement à W‑Sn‑(Cu) de Panasqueira (Portugal). Dissertation, University of Orléans. 2018. Lebocey J. Géologie et gîtologie du district du Pontgibaud (Puy‑ de‑Dôme). In: Bayle LD, De Ascençao Guedes R, Gol D, editors. Le Règne Minéral‑hors série XIX ‑2013. p 35–42. Ledésert B. Fracturation et paléocirculations hydrothermales. Application au granite de Soultz‑sous‑Forêts, Dissertation, Université de Poitiers. 1993. López D‑L, Smith L. Fluid flow in fault zones: analysis of the interplay of convective circulation and topographically driven groundwater flow. Water Resour Res. 1995;31(6):1489–503. Lowell R‑P. Circulation in fractures, hot springs, and convective heat transport on mid‑ ocean ridge crests. Geophys J Int. 1975;40(3):351–65. Lucazeau F. Flux de chaleur, production de chaleur et évolution géodynamique récente du Massif Central Français. Dis‑ sertation, Université des sciences et techniques du Languedoc, Montpellier. 1981. Lucazeau F, Vasseur G, Bayer R. Interpretation of heat flow data in the French Massif Central. Tectonophysics. 1984;103(1):99–119. Lucazeau F, Vasseur G. Heat flow density data from France and surrounding margins. Tectonophysics. 1989;164(2):251–8. Magri F, Möller S, Inbaretal N. 2D and 3D coexisting modes of thermal convection in fractured hydrothermal systems— implications for transboundary flow in the Lower Yarmouk Gorge. Mar Pet Geol. 2016;78:750–8. Malkovsky V‑I, Magri F. Thermal convection of temperature ‑ dependent viscous fluids within three ‑ dimensional faulted geothermal systems: estimation from linear and numerical analyses. Water Resour Res. 2016;52:2855–67. Malkovsky V‑I, Pek A‑A. Timescales for reaching steady‑state fluid flow in systems perturbed by formation of highly permeable faults. Geofluids. 2004;4:253–8. Manning C‑E, Ingebritsen S‑E. Permeability of the continental crust: implications of geothermal data and metamorphic systems. Rev Geophys. 1999;37(1):127–50. Mareschal J‑ C, Jaupart C. Energy budget of the Earth. In: Gupta HK (ed) Encyclopedia of solid earth geophysics. Dordrecht: Springer; 2011. p. 285–91. Martel C, Champallier R, Prouteau G, Pichavant M, Arbaret L, Balcone‑Boissard H, Scaillet B. Trachyte phase rela‑ tions and implication for magma storage conditions in the Chaîne des Puys (French Massif Central). J Petrol. 2013;54(6):1071–107. Duwiquet et al. Geotherm Energy (2019) 7:33 Page 29 of 29 Mazurek M. Geological and hydraulic properties of water‑ conducting features in crystalline rocks. Hydrogeology of crystalline rocks. Dordrecht: Springer; 2000. p. 3–26. Mckenna J‑R, Blackwell D ‑D. Numerical modeling of transient basin and range extensional geothermal systems. Geother ‑ mics. 2004;33(4):457–76. Meredith P‑ G, Main I‑ G, Clint O‑ C, Li L. On the threshold of flow in a tight natural rock. Geophys Res Lett. 2012;39:4. Moeck I‑ J. Catalog of geothermal play types based on geologic controls. Renew Sustain Energy Rev. 2014;37:867–82. Mordensky S‑P, Villeneuve M ‑ C, Kennedy B‑M, Heap M ‑ J, Gravley D‑M, Farquharson J‑I, Reuschlé T. Physical and mechani‑ cal property relationships of a shallow intrusion and volcanic host rock, Pinnacle Ridge, Mt. Ruapehu, New Zealand. J Volcanol Geothermal Res. 2018;359:1–20. Micarelli L, Benedicto A, Wibberley CAJ. Structural evolution and permeability of normal fault zones in highly porous carbonate rocks. J Struct Geol. 2006;28(7):1214–27. Murphy H‑D. Convective instabilities in vertical fractures and faults. J Geophys Res. 1979;84:6121–30. Négroni J‑M. Le district de Pontgibaud cadre Géologique Evolution Structurale et Métallogénique. Dissertation, Univer ‑ sité de Clermont‑Ferrand. 1981. Neuman S‑P. Generalized scaling of permeabilities: validation and effect of support scale. Geophys Res Lett. 1994;21(5):349–52. Oda M. An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses. Water Resour Res. 1986;22(13):1845–56. Person M, Hofstra A, Sweetkind D, Stone W, Cohen D, Gable C‑ W, Banerjee A. Analytical and numerical models of hydro‑ thermal fluid flow at fault intersections. Geofluids. 2012;12(4):312–26. Pribnow D, Schellschmidt R. Thermal tracking of upper crustal fluid flow in the Rhine Graben. Geophys Res Lett. 2000;27(13):1957–60. Rabinowicz M, Boulègue J, Genthon P. Two and three‑ dimensional modeling of hydrothermal convection in the sedi‑ mented middle valley segment, Juan de Fuca ridge. J Geophys Res. 1998;103(B10):24045–65. Ranjram M, Gleeson T, Luijendijk E. Is the permeability of crystalline rock in the shallow crust related to depth, lithology or tectonic setting? Geofluids. 2015;15(1–2):106–19. Roche V, Sternai P, Guillou‑Frottier L, Menant A, Jolivet L, Bouchot V, Gerya T. Emplacement of metamorphic core com‑ plexes and associated geothermal systems controlled by slab dynamics. Earth Planet Sci Lett. 2018a;49:322–33. Rowland J‑ V, Sibson R‑H. Structural controls on hydrothermal flow in a segmented rift system, Taupo Volcanic Zone, New Zealand. Geofluids. 2004;4(4):259–83. Roche V, Bouchot V, Beccaletto L, Jolivet L, Guillou‑Frottier L, Tuduri J, Tokay B. Structural, lithological, and geodynamic controls on geothermal activity in the Menderes geothermal Province ( Western Anatolia, Turkey). Int J Earth Sci. 2018b;104:1–28. Sausse J. Caractérisation et modélisation des écoulements fluides en mi‑ lieu fissuré. Relation avec les altérations hydro ‑ thermales et quantification des paléo ‑ contraintes. Dissertation, Université Henri Poincaré ‑Nancy I. 1998. Sardini P, Ledésert B, Touchard G. Quantification of microscopic porous networks by image analysis and measurements of permeability in the Soultz‑sous‑Forêts granite (Alsace, France). In: Jamtveit B, Yardley BWD, editors. Fluid flow and transport in rocks Mechanisms and effects. New York: Chapman & Hall; 1997. p. 171–89. Smith L, Chapman DS. On the ther ‑ mal effects of groundwater flow: 1. regional scale systems. J Geophys Res. 1983;88(1):593–608. Siebenaller L, Boiron M‑ C, Vanderhaeghe O, Hibsch C, Jessel M, André‑Mayer A‑S, Photiades A. Fluid record of rock exhumation across the brittle‑ ductile transition during formation of a Metamorphic Core Complex (Naxos Island, Cyclades, Greece). J Metamorph Geol. 2013;31:313–38. Sobolev S, Zeyen H, Granet M, Achauer U, Bauer C, Werling F, Altherr R, Fuschs F. Upper mantle temperatures and lithosphere‑asthenosphere system beneath the French Massif Central constrained by seismic, gravity, petrologic and thermal observations. Tectonophysics. 1997;275(1–3):143–64. Sonney R, Vuataz F‑D. Numerical modelling of Alpine deep flow systems: a management and prediction tool for an exploited geothermal reservoir (Lavey‑les‑Bains, Switzerland). Hydrogeol J. 2009;17(3):601–16. Süer S, Güleç N, Mutlu H, Hilton D‑R, Çifter C, Sayin M. Geochemical monitoring of geothermal waters (2002–2004) along the North Anatolian Fault Zone, Turkey: spatial and temporal variations and relationship to seismic activity. Pure Appl Geophys. 2008;165(1):17–43. Taillefer A, Guillou‑Frottier L, Soliva R, Magri F, Lopez S, Courrioux G, Le Goff E. Topographic and faults control of hydro ‑ thermal circulation along dormant faults in an orogen. Geochem Geophys Geosyst. 2018;19:4972–95. Taillefer A, Soliva R, Guillou‑Frottier L, Le Goff E, Martin G, Seranne M. Fault ‑related controls on upward hydrothermal flow; an integrated geological study of the Têt fault, eastern Pyrénées (France). Geofluids. 2017;2017:19. Vasseur G, Gable R, Feuga B, Bienfait G. Groundwater flow and heat flow in an area of mineral springs. Geothermics. 1991;20(3):99–117. Verney N. Importance du contexte géologique dans la signature chimique et isotopique (Carbone 13) des sources hydrothermales. Applications aux émergences du Puy de Dôme. Université Blaise Pascal, MsC thesis. 2005. Violay M, Heap M‑ J, Acosta M, Madonna C. Porosity evolution at the brittle‑ ductile transition in the continental crust: implication for deep hydro‑ geothermal circulation. Sci Rep. 2017;7(1):7705. Vigneresse J‑L, Tikoff B, Améglio L. Modification of the regional stress field by magma intrusion and formation of tabular granitic plutons. Tectonophysics. 1999;302(3–4):203–24. Walter B. Réservoir de socle en contexte extensif: genèse, géométries et circulations de fluids: exemples du rift intraconti‑ nental du lac Albert (Ouganda) et de la marge proximale d’Ifni (Maroc). Dissertation, Université de Lorraine. 2016. Zeyen H, Novak O, Landes M, Prodehl C, Driad L, Hirn A. Refraction‑seismic investigations of the northern Massif Central (France). Tectonophysics. 1997;275(1–3):99–117. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Journal

Geothermal EnergySpringer Journals

Published: Nov 9, 2019

References