THM modeling of hydrothermal circulation at Rittershoffen geothermal site, France

The Rittershoffen deep geothermal project located 6 km east from Soultz-sous-Forts EGS site (France) includes a doublet GRT-1 and GRT-2 to exploit the geothermal resource at the sediments–granite transition where higher temperatures than those of Soultz-sous-Forêts have been measured. Detailed stratigraphic and geophysical data, temperature logs, and tracer surveys have been collected. However, no reservoir model, integrating large-scale geophysical measurements, exists for this site. We developed a reservoir model in two dimensions (10 km × 5 km) based on a finite element method. It includes thermo–hydro–mechanical (THM) coupling and extended brine properties. A representative elementary volume of 100 m is assumed to homogenize the fault network complexity at small scales. A back analysis is performed to obtain large-scale rock properties using GRT-1 temperature log and regional stress-depth profiles. The inverted large-scale properties are consistent with their counterparts measured at the laboratory scale. The bottom of the hydraulic cap rock is 1.2 km ± 0.1 km deep. It is shallower than the discontinuity of the thermal gradient. Hydrothermal convection cells are 2.7 km high which is larger than that previously proposed. A very good fit of the GRT-1 temperature log is obtained using our simplified two-dimensional THM model with four homogenized units at a 100 m scale. The comparison between Rittershoffen and Soultz-sous-Forêts models highlights many similarities in terms of rock properties, decoupling of hydraulic and thermal cap rocks and temperature spatial variability (about 50 °C). Predictions of the relationship between reservoir temperature and surface thermal gradients are proposed for future explorations.


Introduction
In the Upper Rhine Graben (URG), the geothermal gradient is unusually high at the near surface (more than 100 • C/km in the first kilometer in depth, whereas the main value in Europe is 30 • C/km) (Haas and Hoffmann 1929).
The high underground temperatures in the URG make the region the most studied one in Europe for geothermal applications Olasolo et al. 2016;Huenges and Ledru 2011). Geothermal projects, such as the well-known pilot research Soultzsous-Forêts site, are based on enhanced geothermal system (EGS) technology. The EGS concept consists in increasing the reservoir permeability using hydraulic, thermal and/or chemical stimulations and then forcing a circulation of the natural brines in the deep wells taking advantage of the thermal anomaly related to the large-scale hydrothermal system in fractured rocks (Tester et al. 2006;Schindler et al. 2010;Gérard et al. 2006;Schill et al. 2017). After the development (during more 25 years) of the Soultz-sous-Forêts pilot site as an EGS demonstrator, a new industrial project at Rittershoffen was initiated in 2011 and operated in 2016 (Baujard et al. 2015;Genter et al. 2015). The Rittershoffen site is located 6 km east from Soultz-sous-Forêts in Northern Alsace. The project is based on a geothermal doublet, GRT-1 and GRT-2, drilled around 2.6 km deep to intersect the normal Rittershoffen fault and its associated fracture network at the interface between sediments and granite. Structural and stratigraphic studies (Aichholzer et al. 2016;Hehn et al. 2016;Vidal et al. 2016a), temperature logs (Baujard et al. , 2017 as well as seismic (Gaucher et al. 2013;Maurer et al. 2015;Lengliné et al. 2017) and geochemical surveys (Dezayes et al. 2014; have been already established to identify permeable zones and the hydraulic connections between GRT-1 and GRT-2 wells. The knowledge of the Rittershoffen site is also completed by the huge database collected in the nearby Soultz-sous-Forêts site. Five permeable zones have been currently identified in the granite reservoir, but none in the sediments (Vidal et al. 2016b).
As already demonstrated for the Soultz-sous-Forêts research site, numerical modelings can provide significant insights to better understand the hydrothermal circulation or the rock physics of EGS reservoirs (Jain et al. 2015;Kolditz and Clauser 1998;Pruess 1990;Sanyal 2000;Tomac 2017). Numerical reservoir models can typically be classified according to different aspects (Willis-Richards and Wallroth 1995): (i) the description of the complex fracture network geometries they integrate via stochastic distribution (Baujard and Bruel 2006;Cacas et al. 1990) or regular grids (Watanabe and Takahashi 1995;Kohl and Mégel 2007;Willis-Richards et al. 1996), (ii) the analysis of detailed physical processes such as full thermo-hydro-mechanical-chemical (THMC) couplings they account for (Kohl et al. 1995;Gelet et al. 2012;Diersch and Kolditz 1998;Bachler and Kohl 2005;McDermott et al. 2006a, b). Thermo-hydro-mechanical (THM) modeling based on a homogenized description of the reservoir has been recently presented in Magnenet et al. (2014) for the Soultz-sous-Forêts site to describe the natural hydrothermal circulation.
In the continuity of this work, the current study aims at proposing a new model of the large-scale hydrothermal circulation in the recent Rittershoffen EGS site. The numerical modeling is based on the current local geological and geophysical knowledge from Rittershoffen, but also the closeby Soultz-sous-Forêts site. The 2D reservoir model includes all major THM couplings. The equations governing THM processes are solved by a finite element approach using the Code_Aster software. The reservoir is homogenized at the scale of 100 m. The fluid rheology (e.g., density, viscosity, heat capacity) is considered as dependent on temperature and fluid pressure as shown by laboratory measurements (Zaytsev and Aseyev 1992;Kestin et al. 1981;Rowe and Chou 1970). We include different a priori settings: (i) the main geological structures of the sedimentary cover (Aichholzer et al. 2016) and the basement (Vidal et al. 2016a), (ii) the temperature-depth profiles through the deep boreholes GRT-1 and GRT-2 (Baujard et al. 2017), (iii) the distribution of the natural radioactivity (Rummel 1992;Pribnow et al. 1999;Pribnow and Schellschmidt 2000), (iv) the regional stress-state (Evans et al. 2009;Cornet et al. 2007;Valley 2007), (v) the rock properties and their up-scaling (GeORG13), and (vi) the geochemical data obtained from brine samples . We proceed to a back analysis to find the reservoir parameters from the temperature and stress logs. We finally compare the results with the insights from the Soultz-sous-Forêts site and discuss the impact of the geological settings and large-scale fault on the hydrothermal circulation, the location of the hydraulic and thermal cap rocks, the lateral variability of the reservoir temperature and the link between the geothermal gradient and the reservoir temperature at a depth of 2.0 km. a b Fig. 1 a 2D conceptual model of the geology at Rittershoffen. The sedimentary cover is investigated from geological studies from Vidal et al. (2017) and Aichholzer et al. (2016). The granite below 2.5 km deep is assumed to be the same as in the Soultz-sous-Forêts site (Dezayes et al. 2005a, b). b Temperature-depth profiles obtained from logs run in the Rittershoffen wells (GRT-1, GRT-2) after drilling operation over (Baujard et al. 2017). The background colors correspond to the four layers homogenized at the scale of about 100 m and considered in the model ( e 1 , e 2 , e 3 and e 4 correspond to the layer thicknesses, e 1 is inverted during the back analysis; e 2 = 2.2 km -e 1 ; e 1 + e 2 + e 3 = 3.9 km ; e 1 + e 2 + e 3 + e 4 = 5.4 km) Figure 1a shows a representation of the main geological units of the Rittershoffen site. Its geology is similar to those of Soultz-sous-Forêts site except for the thicknesses of the units (Aichholzer et al. 2016). At Rittershoffen, the first 2200 meters consists of a sedimentary cover overlying a granitic basement, whereas the interface sedimentsgranite is at 1400 m for the Soultz-sous-Forêts site. The sedimentary sequence begins by a sandy-clayey Pliocene-Quaternary layer. This layer is thicker by about 530 m at Rittershoffen than at Soultz-sous-Forêts due to a higher erosion rate in the latter. The Pliocene-Quaternary sequence directly lays on clays and marls of the Oligocene age. Contrary to Soultz-sous-Forêts, the Rittershoffen site exhibits a full Grey Series complex and an upper part of the Pechelbronn layers. Below, Eocene formations are composed of two units: ferruginous marls of the red layer and clayey marlstones interbedded with domolite layers (i.e., the dolomitic zone). The dolomitic zone is also thicker in Rittershoffen than in Soultz-sous-Forêts. The Tertiary unit is overlying Jurassic formations that include Dogger black silty clay shales and Lias gray clay formations. The Triassic sequence is the deepest part of the sedimentary cover. The Triassic units are thicker in Rittershoffen than in Soultz-sous-Forêts. It corresponds to Keuper marls and clays, Muschelkalk limestones and marly calcareous dolomites. The last Triassic layer corresponds to the Buntsandstein sandstone. Core studies show that the top of the granitic basement is divided into three parts . From the top to the bottom, it is composed of reddish oxidized granite due to paleo-weathering, hydrothermally altered granite and fine-grained low altered granite. Here, we assume that after 2.5 km in depth, the granitic basement is the same in Rittershoffen as in Soultz-sous-Forêts. The basement until 3.9 km deep is composed of a porphyritic monzo-granite with K-feldspar megacrysts in Soultz-sous-Forêts (Dezayes et al. 2010(Dezayes et al. , 2005c. Located below this is the first transition at about 3.9 km to a biotite and amphibole enriched granite and the second transition at about 4.6 km to a rather different leuco-granite with very finegrained micas.

Geological settings
Ultra borehole images (UBIs) logs and geochemical analyses have been performed to investigate the structural properties of the fracture networks in Rittershoffen (Vidal et al. 2016a;Dezayes et al. 2014;Vidal et al. 2017). As in Soultz-sous-Forêts, two main natural fracture systems have been identified (Dezayes et al. 2014). The first is composed of closely connected meso-fractures. The second is a set of large fractures crossing the former system. From structural analysis, the main set of fractures is oriented around N15-20 • E with a dip of 80 • W in GRT-1, but more scattered in GRT-2. In the sediments, the density of fractures is about 0.33 fractures per meter. Among eight major fracture areas (i.e., with a thickness higher than 1 cm), only one fracture cluster is considered as permeable (Vidal et al. 2016b). It has an orientation of N20 • E with a dip of 85 • W. In contrast, the top of the basement is highly fractured with about 2.51 fractures per meter, even more than in Soultz-sous-Forêts (with 0.65 fractures per meter). Four permeable fracture areas have been observed among 11 major fracture zones. The fracture zones in the granite have the same orientation as in the Soultz-sous-Forêts site, parallel to the regional orientation. Vallier et al. Geotherm Energy (2018) 6:22 Temperature profiles and hydrothermal circulation Figure 1b illustrates the temperature-depth profiles in GRT-1 and GRT-2 (Baujard et al. 2017). From the surface to the top of the Muschelkalk, the geothermal gradient is constant (around 85 • C km −1 ) in both wells. The value is slightly lower than in Soultz-sous-Forêts where it is about 110 • C km −1 . Below, the geothermal gradient suddenly declines about 30 times, i.e., around 3 • C km −1 in GRT-1 and about 18 • C km −1 in GRT-2 at the time of the measurement. The difference in geothermal gradient between the two wells is explained by the thermal non-equilibrium of the GRT-2 well . Some local temperature perturbations have been recorded in the profiles between 1500 and 2700 m in depth. They are commonly considered as evidence of hydrothermal circulation through fracture zones in particular around 1650 and 2350 m deep in GRT-1. The temperature evolution with depth is unknown below the bottom of the wells. Temperature logs in Soultz-sous-Forêts suggest that the geothermal gradient tends to the average Central European gradient (around 30 • C km −1 ) in the deep granitic basement Pribnow et al. 1999). From the surface to the top of the Muschelkalk, the linear temperature trend suggests that the thermal state is purely diffusive. The sharpness of the transition at 1.65 km deep would imply that the upper sediments behave as a hydraulic cap rock. Below, the very low geothermal gradient is classically interpreted as the result of an active and large-scale hydrothermal convection within the fractured granite. Indeed, an intense circulation of the native brine has been evidenced in the reservoir from geochemical analyses of the in situ fluid in GRT-1 and GRT-2 (Dezayes et al. 2014;. The low tracer concentrations in GRT-2 during an injection test in GRT-1 highlight a weaker connection between the wells than in Soultz-sous-Forêts . However, we will show in the present study that this interpretation is oversimplified and has to be reconsidered. The fluid analysis indicates that natural brine has the same salinity, pH and chemical characteristics as in Soultz-sous-Forêts. The contribution of natural granite radioactivity to the origin of the thermal anomaly in the Upper Rhine Graben has been studied from core analyses at the Soultzsous-Forêts site (Rummel 1992;Pribnow et al. 1999;Pribnow and Schellschmidt 2000;Baillieux et al. 2013). The production rates are typically of the order of 0.1, 1.0, 5.0 µ W m −3 , respectively, for the upper sediments, the Buntsandstein sandstone and the two granites (Kohl 2000). We will assume that they are similar at Rittershoffen.

Rock physics
No direct laboratory measurements of the sediments and granite properties are available for the Rittershoffen site. This is why the rock properties at Rittershoffen are assumed to be the same as at the vicinity of the Soultz-sous-Forêts site. Table 1 presents a synthetic review of the relevant rock properties obtained either from laboratory measurements on core samples or from geophysical investigations. Thermal conductivities vary between 1.1 and 5.9 W m −1 K −1 (GeORG 2013) and permeabilities between 1.0 × 10 −20 and 3.2 × 10 −12 m 2 with a large variability which will be used as prior distribution for the back analysis (Hettkamp et al. 1999;Kohl 2000;Bar 2012;GeORG 2013;Magnenet et al. 2014;Griffiths et al. 2016;Heap et al. 2017).

Toward a large-scale reservoir model at Rittershoffen
The goal of the study is to build the simplest THM numerical model that is consistent with the main characteristics of the Rittershoffen site. The model does not aim to describe all the complexity of the geology or the deterministic details of the faults networks ( Fig. 1a). As sketched in Fig. 1b, the whole sedimentary cover is split into two horizontal homogenized units: the upper sediments and the lower sediments. The depth of the transition between the two units is taken as a parameter to be adjusted (named e 1 ) during the parameter back analysis. The basement is also split into two units: the upper granites and the lower granites. Due to the lack of direct knowledge on the deep granitic basement in Rittershoffen, the transition between the two units is set at a depth e 1 + e 2 + e 3 = 3.9 km as in Soultz-sous-Forêts site. The transition between the sediments and the granite is also set at e 1 + e 2 = 2.2 km in depth. We assume that the radiogenic sources are set at 0.1, 1.0, 5.0 µW m −3 , respectively, for the upper sediments, the lower sediments and the two granites (Kohl 2000).

Governing equations of the THM model
All units are homogenized as a porous medium, fully saturated with a single-phase brine and including thermo-hydro-mechanical (called THM) coupling as developed in Coussy (2004). Here are our main assumptions: • A small perturbation assumption is made and solid grains are considered to remain in the thermo-elastic regime; • The Cauchy stress tensor σ is split into two contributions: an effective stress σ ′ and a hydraulic stress σ p 1 ( 1 being the unit tensor); Total specific mass, r 0 (kg m −3 ) Young's modulus, E Heat source produc- • The thermodynamic flows (effective Cauchy stress σ ′ , water surface mass flow M w , heat flow q ) are lineary related to thermodynamic forces (linearized strain ǫ , gradient of pore pressure ∇p w , gradient of temperature ∇T ), but with coefficients that may depend on temperature, porosity (denoted here φ ) or pore pressure. Hence, most of the homogenized properties (such as the specific heat at constant stress and the thermal conductivity) appearing in Hooke's law, Darcy's law, and Fourier's law of the porous materials are considered as functions of the form f (φ, p w , T ) by using classical mixing laws; • Following the approach of Magnenet et al. (2014) and Vallier et al. (2016), we consider a rheology of brine that is extrapolated from experimental results for artificial brines at different salinities (NaCl) (Zaytsev and Aseyev 1992;Kestin et al. 1981;Rowe and Chou 1970). More specifically, we assume that the natural brine is equivalent to a pure NaCl solution with a mean specific mass content of 100 g L −1 . The retained mathematical expressions of the brine properties (density, dynamic viscosity, thermal dilatation, thermal conductivity) are given in Table 2.
The whole set of notations as well as a detailed presentation of the governing equations is presented in Appendix 1.

The finite-element model
The governing equations are solved by using the open-source finite element solver Code_ Aster, (EDF, 2016;Zienkiewicz et al. 2013) in which specific developments were added to account for the heat sources induced by the radioactivity of rocks, the nonlinear rheology of brine and the search of the stationary solutions. Figure 2 shows the considered cross section of the four idealized units from our Rittershoffen reservoir model (Fig. 1b).

Parameter Expression Coefficients
This idealized cross section is 5.35 km in height and 10 km in width and contains no large-scale fault. Below, we denote by x the horizontal direction, z the out-of-plane horizontal direction and y the vertical one. The typical size of elements is 100 m × 100 m. Consequently, the rheological properties taken as inputs as the model are homogenized at this length scale. In Fig. 2, the boundary conditions used in our numerical simulations are also summarized: • Temperatures are, respectively, maintained at 10.0 and 213.0 • C on the upper and lower boundaries. The lateral boundaries are taken as adiabatic. • A fluid pressure of 0.1 MPa (i.e., the value of atmospheric pressure) is imposed on the upper boundary. The other boundaries are assumed to be impermeable. • The normal displacement is nil on the lower and lateral boundaries. The upper boundary is stress free.
A Euler implicit scheme is used for the time integration of nodal mechanical displacements, pore pressures and temperatures. The nonlinear system is solved by the Newton-Raphson method. For the initial conditions, a constant and uniform temperature of 10.0 • C is assumed. The fluid pressure field is also assumed to be constant at 0.1 MPa. To ensure the convergence of the process, the computation has been divided into three steps (Magnenet et al. 2014): (i) during a short time period of 1000 years, the boundary conditions and gravity are progressively applied ; (ii) next, during 100,000 years, the system freely evolves along constant boundary conditions ; (iii) in one last increment, the system reaches a steady state by cancelling the nonstationary terms from the constitutive equations.

Inverse method
In this study, some geometrical and rheological parameters are estimated by back analysis (see next section). To do it, Code_Aster has been coupled to the Parameter ESTimation (PEST) software (Doherty 2005).
Numerous inverse methods have been established to carry out back analysis of the rock properties: Monte-Carlo methods, Bayesian approaches which associate probability distribution for each parameter (Vogt 2012;Kosack et al. 2011;Tarantola 2004), the neighborhood algorithm based on random generation of new parameters (Sambridge 1999) or the genetic algorithm (Pérez-Flores and Schultz 2002).
Here, the PEST back analysis software is an implementation of the so-called Levenberg-Marquardt algorithm which minimizes an "error function"-typically the L 2 -norm of the difference between model and observations-with respect to a chosen set of parameters. Each parameter p is taken from a uniform a priori distribution called prior distribution in the range [p min , p max ] chosen to be wider than experimental values (see Table 1).
The main benefit of our back analysis is that the inversion procedure based on deterministic method is less numerically intensive than the stochastic methods. However, the inversion procedure is sensitive to the initial conditions of the back analysis such as the prior distributions of the rock properties. Nonetheless, the prior distributions are well constrained thanks to databases from Rittershoffen and Soultz-sous-Forêts sites.
In the current study, we aim to estimate four rock properties : permeability, thermal conductivity and elastic moduli (i.e., Young's modulus and Poisson's ratio) by inverting the observed temperature and stress profiles with depth using the THM model as the forward model. During this back analysis, the four rock properties are explored in the prior uniform distributions described in Table 3. The prior distributions, wider than the variability of the data, allow the inversion to explore all the relevant experimental values for the Rittershoffen site. Other rock properties described in Table 1 are set as constant during the back analysis. Their values are summarized in Table 4.    During the PEST inversion, several assumptions have been made: • The depth of the granitic basement is taken as the same as in Soultz-sous-Forêts.
More precisely, the interface between the upper and lower granites is set at the depth e 1 + e 2 + e 3 = 3.9 km. We also set the interface between sediments and granite at e 1 + e 2 = 2.2 km. • The ranges of values of rock properties are assumed to be the same as in Soultz-sous-Forêts (see Table 1). • The observed temperature profile is assumed to be at the location of the surface maximum heat flux. This assumption is consistent with surface temperature maps in the Rittershoffen area (Haas and Hoffmann 1929). In practice, the numerical temperature profile has then been calculated at the side of a convective cell where the Darcy's velocity is purely ascending. The stress-depth profiles are taken from the same position.

Results of the back analysis
The goal of the study is to better understand the physics of the Rittershoffen reservoir. To address this issue, a back analysis confronting our THM model with the observed temperature and stress depth profiles is carried out. To reproduce the GRT-1 a b c temperature-depth profile, nine parameters are estimated: the thickness of the first geological layer e 1 , and the permeabilities K int,i and thermal conductivities i of the four layers i = 1 . . . 4. An excellent fit was found between the simulated temperature-depth profile and the measured GRT-1 one. The main trend of the observed T-log (except for the sharpness observed around 1.65 km deep), the temperature gradients and the reservoir temperature around 160 • C are well reproduced. The temperature and Darcy's velocity maps of the best model are displayed in Fig. 3. The convection cells have a width of about 3.0 km and a height of 2.7 km. The maximum of Darcy's velocity is about 20.0 cm year −1 . The hydraulic cap rock corresponds to a contrast of permeability and isolates the hydrothermal circulation as convection below. The depth of the top of the convection cells is the depth of the hydraulic cap rock. By analogy, the thermal cap rock corresponds to a contrast of thermal conductivity. Here, the hydraulic cap rock (i.e., the top of the convection cells) appears shallower than expected from previous interpretations of the experimental temperature-depth profiles from GRT-1 and GRT-2 wells (Baujard et al. 2017). Previous studies suggest that the hydraulic cap rock should correspond to the top of the Muschelkalk formations at a depth of 1.65 km. However, the hydraulic cap rock is obtained in the best model at a depth of 1.2 km ± 0.1 km and is associated to a very good fit of the observed GRT-1 T-log. As illustrated in Fig. 3, the depth is close to the depth for the bottom of the Tertiary formations and not the depth of the transition between the sediments and the granitic basement. Interestingly, the bottom of the cap rock does not correspond to the breaking point of GRT-1 log located at the Muschelkalk top. Darcy's velocities in our model are of the same order of magnitude as the estimated values from previous numerical studies in the Upper Rhine Graben (Clauser 1990;Kohl 2000;Guillou-Frottier et al. 2013) and hydraulic tests (Bari et al. 1998). However, the Darcy's velocity is slightly higher than the one predicted in the Buntsandstein sandstone (between 1.8 and 2.2 km in depth) of the order of 5-10 cm year −1 by (Guillou-Frottier et al. 2013). This may be due to the high large-scale permeability found in the lower sediments and upper granites, higher than the one found for a similar model of Soultz-sous-Forêts (see Fig. 4). Figure 3 shows the variation of the vertical component of the surface heat flow along the x-axis in the model. The average surface heat flow is 128 mW m −2 ± 12 mW m −2 . This value is in accordance with that found in the Upper Rhine Graben Region, 150 mW m −2 in Soultz-sous-Forêts (Bächler et al. 2003;Pribnow and Schellschmidt 2000;Clauser et al. 2002).
The range of laboratory values (as hatched zones), the prior distributions (as dashed lines) and the inverted properties (as thicked lines) are plotted in Fig. 4. The inverted permeabilities are, respectively, around 1.0 × 10 −16 m 2 and 8.0 × 10 −16 m 2 in the upper sediments (above 1.2 km of depth) and lower granites. The permeability is 1.6 × 10 −14 m 2 in the lower sediments (below 1.2 km deep) and upper granites. The permeabilities are in good agreement with laboratory values. Concerning the thermal conductivities, they are, respectively, 1.4 and 3.1 W m −1 K −1 in the sediments and granites, which is also consistent with laboratory values.
Two cap rocks can be identified from the estimated vertical profiles of K int and . The hydraulical cap rock is associated at its base to the high contrast of permeability and the thermal cap rock to the discrepancy of thermal conductivity. Here, the bottom of the hydraulical cap rock (i.e., the top of the convection cells, see Fig. 3) is identified at the interface between the upper and lower sediments. The contrast of permeability is associated here with the high fracture density in the lower sediments and the granite compared to the upper sediments (Vidal et al. 2016b. The change of the rock property does not correspond to a contrast in terms of lithology. On the contrary, the contrast of thermal conductivity is located at the interface between the sediments and the granitic basement. The whole sedimentary cover associated with a lower thermal conductivity than in the basement contributes to a thermal blanketing of the insulating sediments. The effect has been already identified as a key factor to explain the higher geothermal gradient at depth than the average European one (Freymark et al. 2017;Scheck-Wenderoth et al. 2014). The discrepancy in terms of depths between the permeability and thermal conductivity contrast highlights a decoupling of the cap rocks: the whole sedimentary cover corresponds to the thermal cap rock, whereas only the upper sediments behave as a hydraulic cap rock.
In our computation, the regional Rittershoffen stress state is assumed to be similar to that of the Soultz site as proposed in Baujard et al. (2017). The observed trends of principal stress values with depth, established for the Soultz site (Evans et al. 2009), are then used for the Rittershoffen reservoir. They are presented in Appendix 2 as Eqs. 19, 20 and 21. To reproduce these observed trends of the stress state, a back analysis of the elastic

Fig. 4 Comparison of permeability estimated by back analysis (left) and thermal conductivity (right) between
Rittershoffen and Soultz-sous-Forêts sites. The shadow zones correspond to the range of experimental values (permeabilities in green, and thermal conductivity in red) (see Table 1). The dashed lines correspond to the prior distributions for the back analysis (see Table 3). Background colors correspond to geological layers moduli (Young's modulus and Poisson's ratio) is performed. The model is assumed to be oriented along the direction of the maximum horizontal principal stress. Figure 5 presents the best fits by the simulated principal stress-depth profiles of the observed trends. The inverted Young's modulus and Poisson's ratio-depth profiles (as thick lines) are plotted in the new Fig. 6. The range of laboratory values are plotted as hatched zones, and the prior distributions as dashed lines. The inverted Young's moduli are respectively 15 and 25 GPa for the sediments and for the granites. The inverted Young's moduli are small compared to the laboratory values in particular for the granites. They are the longterm moduli, different from the short-term measurements at the laboratory scale. The inverted Poisson's ratio are, respectively, 0.23, 0.25 and 0.20 for the sediments, the upper and the lower granites, in accordance with experimental values.

Influence of the large-scale fault in the THM model
In the present study, the hydrothermal circulation is assumed to be mainly driven by fractures and faults networks having a dimension less than the representative elementary volume (REV) with a size of 100 m. The hypothesis is different from many Soultzsous-Forêts modeling studies, which consider that large-scale faults (i.e., larger than 100 m) contribute dominantly to the thermal state and accordingly have to be included in the reservoir model (Baujard and Bruel 2006;Kohl and Mégel 2007;Kohl 2000).
To support our assumption, a large-scale fault has been included in our reservoir model. It corresponds to the Rittershoffen fault, a major fault zone with a north-south strike. It extends from the surface to 3.5 km deep (GeORG 2013), with a thickness of  (Baujard et al. 2017; ii) 74 • in a regional scale best fitting plane of the seismic cloud from induced seismicity (Lengliné et al. 2017;iii) 83 • from small-scale acoustic image logs (Vidal et al. 2016a). The difference between the estimates may be linked to the scale discrepancy between the three approaches: the structural model is built at the 5 km scale, whereas the acoustic logs are performed at the scale of 0.2 m. Simulations have been carried out with the different estimated dips: 45 • , 74 • and 83 • to the west. Figure 7 shows the comparison of the temperature-depth profiles between the cases with and without fault and Fig. 8 shows the stress-depth profiles. The maximum of temperature difference between the cases with and without the Rittershoffen fault is about 6 • C at 2.0 km deep. Temperatures are slightly higher in the sediments and weaker in the granite after adding the Rittershoffen fault. Differences in the fault dip have also negligible impact on the temperature distribution (less than 2 • C at 2.0 km in depth). Concerning the stress-depth profiles, the difference is even less noticeable (about 1 MPa) for the maximum horizontal stress at the transition between the upper and lower granites. All the stress-depth profiles still fit the experimental profiles of the three principal stress components with depth. Figure 9 shows the distribution of the Darcy's velocities in the Fig. 6 Comparison of Young's modulus (left) and Poisson's ratio (right) between Rittershoffen and Soultz-sous-Forêts sites. The shadow zones correspond to the range of experimental values (Young's modulus in yellow, and Poisson's ratio in gray) (see Table 1). The dashed lines correspond to boundaries of the prior distributions for the back analysis (see Table 3). Background colors correspond to geological units model including the fault with a dip of 45 • . A circulation occurs upward along the fault with a maximum of Darcy's velocities of about 26.0 cm year −1 , slightly higher than in the case without fault. However, the system of convection cells stays undisturbed; their numbers and sizes are the same after inclusion of the fault. In conclusion, the simulations adding a large-scale fault in the model do not show an important influence on the simulated temperature and stress-depth profiles justifying the main assumption of our approach.

Influence of the different couplings in the THM model
The claim of a full thermo-hydro-mechanical analysis stands in our modeling approach by the inclusion of several coupling between the three physical processes. The material properties depend on field variables such as pressure, temperature and porosity (see Table 2). Concerning the poroelastic behavior of the system, the Cauchy stress is assumed to be split into two components with one being the hydraulic stress with its variation depending on one of the fluid pressures (see Eq. 8 in Appendix 1). The poroelasticity is also described by the incremental variation of porosity depending on the To study the effect of the different couplings on the thermal state, simulations have been carried out after canceling them. Figure 10 shows the temperature-depth profiles and the associated maps of temperatures for these simulations. We can observe that the cancellation of the dependence of the thermics on the mechanical and hydraulic processes (HM) or of the hydraulics on the mechanics and thermics (TM) leads both to a large-scale diffusive case. The deactivation of the mechanical effect on the thermal and hydraulic processes (TH) highlights a convection system, but different from the full THM model. It includes only two convective cells larger than the ones from the THM model and the reservoir temperature is slightly higher than in the previous model. To conclude, cancelling the mechanical or hydraulic parts of the THM model does have a clear impact on the thermal regime and the hydrothermal circulation.

Comparison with the back analysis for the closeby Soultz-sous-Forêts site
A two-dimensional THM model has been also developed for the Soultz-sous-Forêts site (Vallier et al. 2017). Permeabilities, thermal conductivities, Young's moduli and Poisson's ratios have been estimated to reproduce the GPK-2 temperature log and the total stress profiles via a similar back analysis. In Soultz, the best model highlights a hydrothermal circulation below a shallow bottom of the hydraulic cap rock at a depth of 100 m, whereas it is at a depth of 1200 m in Rittershoffen. Figure 4 illustrates the estimated In Soultz, the inverted permeabilities are, respectively, 1.0 × 10 −17 m 2 , 3.5 × 10 −15 m 2 , 6.0 × 10 −15 m 2 and 2.5 × 10 −16 m 2 for the upper sediments, the lower sediments (below 100 m deep), the upper and the lower granites. The permeability increases with depth from the upper sediments to the lower sediments and upper granites. Importantly, the permeability of the lower sediments and the upper granite is very similar suggesting that the lithological transition between the sedimentary cover and the granitic basement is not significant for the hydraulic properties. The permeability decreases after 3.9 km in depth.
For the thermal conductivity in Soultz, the values are 3.1 Wm −1 K −1 in the granites and 2.1 Wm −1 K −1 in the sediments. The thermal conductivity increases at the interface between sediments and granitic basement. Then, its value remains constant in the whole granitic basement. The thermal property is controlled by the interface between the sediments and granites on the contrary to permeability.
Profiles for the two properties highlight that both geothermal sites share noticeable similarities. Inverted properties of the best models show common general trends and in particular, the decoupling of behaviors between the thermal and hydraulic cap rocks. Both sites show a blanketing effect from the whole sedimentary cover, whereas the top of the convective cells (i.e., the bottom of the hydraulic cap rock) is at the Fig. 9 Maps of calculated temperatures (background colors) and Darcy's velocities (black arrows) in the model including the Rittershoffen fault with dip of 45 • transition between the upper and lower sediments. The other similarity between both geothermal sites concerns the rock properties. The reservoir permeability is similar between both sites, about 1.0 × 10 −14 m 2 in the granitic basement and the lower sediments. The permeability is also very close in the deep granites, about 6.0 × 10 −15 m 2 . The slightly higher permeability for the Rittershoffen reservoir can be explained by the more important fracture density in the granite compared to Soultz-sous-Forêts (Vidal et al. 2016b).
There is a noticeable discrepancy between the two sites: the thicknesses of the hydraulic cap rocks. The difference is more than 1 km and leads to a discrepancy between the permeabilities in the upper sediments. This difference can be explained by a higher fracture density for the sediments in Soultz-sous-Forêts than in Rittershoffen. This is consistent with the recent stratigraphic studies comparing the two geothermal sites (Aichholzer et al. 2016). Indeed, a more intense fault network has been observed for the sediments in Soultz-sous-Forêts (Aichholzer et al. 2016;Vidal et al. 2016a). Figure 6 also provides a comparison of the estimated moduli for Rittershoffen and Soultz sites. The Young's modulus and Poisson's ratio are very similar in particular for the granitic basement between sites.

Fig. 10
Left: comparison of temperature-depth profiles for four studied cases: in green, the full coupling between thermic, hydraulic and mechanic processes is taken into account; in red, the coupling is only between hydraulics and mechanics; in blue, between thermics and mechanics; and in orange, between thermics and hydraulics. Right: associated maps of temperatures

Comparison with the hydrothermal characterization of the GRT-1 and GRT-2 wells
The study of (Baujard et al. 2017) analyzes the database obtained after stimulation and circulation testing of GRT-1 and GRT-2 wells in the Rittershoffen reservoir. In particular, a double porosity model area has been carried out with AQTESOLV software ) from the hydraulic data of the pumping tests. The reservoir thickness is assumed to be about 500 m and the model includes a fracture area 40 m thick. We aim to compare the results from our modeling approach with their interpretations of production and circulation tests.
Fracture and matrix permeabilities of, respectively, 5.34 × 10 −14 m 2 and 9.2 × 10 −15 m 2 have been evaluated from interpretations of GRT-2 hydraulic testing (Baujard et al. 2017). In our model, the computed permeability for the convective area (i.e., the lower sediments and the upper granites) is about 1.6 × 10 −14 m 2 . This value is between the matrix and fracture permeabilities proposed by Baujard et al. (2017). Our inverted permeability is then consistent with the values from the interpretations of hydraulic tests. Baujard et al. (2017) proposed that the depth of the top of the Muschelkalk formations, i.e., 1.65 km, is the depth of the hydraulic cap rock owing to the transition in the temperature log. In our study, the depth of the hydraulic cap rock has been evaluated at 1.2 km ± 0.1 km after our back analysis, which is significantly above the important change of the temperature gradient. From the interpretation of temperature logs, the height of the convection cells is 1350 m (Baujard et al. 2017). In our model, the simulated convection cells have a height of 2.7 km. They extend shallower into the sediments and much deeper into the granitic basement until 3.9 km in depth assuming that the deep granites in Rittershoffen are similar to the ones from the Soultz-sous-Forêts site. Baujard et al. (2017) has also evaluated the Rayleigh number at Rittershoffen. This calculation aims to confirm that a hydrothermal convection occurs inside the Rittershoffen reservoir. The Rayleigh number is found to be included between 11.1 and 535.7. By using the same definition of the Rayleigh number (Desaive 2002;Turcotte 2014), we found a value of 50. Interestingly our value is included in the range found by Baujard et al. (2017). Knowing that the critical Rayleigh number is about 39.5 (Turcotte 2014), it confirms that a spontaneous convection inside the fractured granite is expected in the Rittershoffen model (Murphy 1979).

Temperature lateral variability
In the prospect of future geothermal exploitation, a precise assessment of the reservoir temperature at depth is required from measurements acquired on the near surface. To address the issue, we aim to bring some insights concerning the link between the lateral variability of the reservoir temperature at 2000 m depth and the one of the geothermal gradient obtained at the near surface. Figure 11 shows different temperature-depth profiles taken at several horizontal positions (every kilometer). Measured T-logs GRT-1 and GPK-2 have been added for comparison. The lateral variability of the reservoir temperature at 2000 m deep is about 40-50 • C. To be noted, the variability of the temperature-depth profiles in the Rittershoffen large-scale model is not enough to reproduce the GPK-2 temperature-depth profile observed at the Soultz-sous-Forêts site even if both geothermal sites share similarities in terms of rock properties.
To better understand the link between the lateral variabilities of the reservoir temperature and the geothermal gradient, Fig. 11 illustrates also the variation of geothermal gradient at the near surface (for the first 200 m depth) along the x-axis. The thermal gradient varies between 76 and 91 • C km −1 . The same periodicity of 6 km is observed for both the geothermal gradient and the reservoir temperature. Interestingly, this periodicity of 6 km can be compared to the distance between Soultz-sous-Forêts and Rittershoffen sites (around 6.5 km).

Conclusion
By using a back analysis confronting a THM model to the temperature and stress profiles observed at Rittershoffen, an excellent fit of the T-log has been found as well for the regional stress-depth trends. The bottom of the hydraulic cap rock (i.e., the top of the convection cells where a contrast of permeability is obtained) is identified at a depth of 1.2 km ± 0.1 km. This depth is close to the bottom of the Tertiary formations and does not correspond to the discontinuity of temperature-depth profile observed in the GRT-1 T-log at 1.65 km deep. The computed permeability is 1.6 × 10 −14 m 2 for the Fig. 11 a Temperature map showing the extracted temperature-depth profiles to study the lateral temperature variability. b Simulated vertical temperature profiles (as red lines) obtained every kilometer along a half-convection cell compared with the observed GRT-1 log from Rittershoffen (as blue dashed line) and GPK-2 log from Soultz-sous-Forêts (as green dashed line). c Thermal gradient obtained from the surface obtained every kilometer along a half-convection cell lower sediments and the upper granite. This highlights that the lithological transition between the sediments and the granitic basement has little influence on the hydraulic property. Contrary to the hydraulic cap rock, the bottom of the thermal cap rock (i.e., the zone of thermal conductivity contrast) occurs at the interface between the sediments and the granite. The thermal conductivity is, respectively, 1.4 and 3.1 W m −1 K −1 in the sediments and granites. This means that the whole sedimentary cover contributes to a blanketing effect and that the thermal and hydraulical cap rocks are decoupled. The large-scale permeabilities, thermal conductivities and elastic moduli are mostly consistent with the values observed at the laboratory scale. They have been compared to the ones obtained from a similar back analysis at the closeby Soultz geothermal site. The permeability, thermal conductivity, Young's modulus and Poisson's ratio have the same general trend with depth and similar values at Soultz and Rittershoffen. Both sites highlight the same decoupling of the hydraulic and thermal cap rocks. The lateral variability of the reservoir temperature at 2.0 km deep is similar between Rittershoffen and Soultzsous-Forêts around 40-50 • C. The same lateral periodicity of 6 km has been found for the geothermal gradient obtained from the near surface and the reservoir temperature. This might lead potentially to a promising tool to assess future geothermal resources. Moreover, further works are currently being done to investigate the potential influence of major faults in the model (e.g., here the Rittershoffen fault). This will allow us to better understand their influence on both thermal and mechanical behaviors inside the reservoir.