Hydro-Thermal Modeling for Geothermal Energy Extraction from Soultz-sous-Forêts, France

The deep geothermal energy project at Soultz-sous-Forêts is located in the Upper Rhine Graben, France. As part of the Multidisciplinary and multi-contact demonstration of EGS exploration and Exploitation Techniques and potentials (MEET) project, this study aimed to evaluate the possibility of extracting higher amounts of energy from the existing industrial infrastructure. To achieve this objective, the effect of reinjecting fluid at lower temperature than the current fluid injection temperature of 70 °C was modeled and the drop in the production wellhead temperature for 100 years of operation was quantified. Two injection-production rate scenarios were considered and compared for their effect on overall production wellhead temperature. For each scenario, reinjection temperatures of 40, 50, and 60 °C were chosen and compared with the 70 °C injection case. For the lower production rate scenario, the results show that the production wellhead temperature is approximately 1–1.5 °C higher than for the higher production rate scenario after 100 years of operation. In conclusion, no significant thermal breakthrough was observed with the applied flow rates and lowered injection temperatures even after 100 years of operation.


Introduction
Geothermal energy is a clean, renewable and low-cost solution for heating and power generation. One of the most challenging problems that humanity is facing is how to mitigate climate change and the anthropogenic emission of carbon dioxide, in order to achieve the target of the Paris agreement, which limits the atmospheric temperature rise to 2 • C or less [1]. Carbon geosequestration is the most desirable solution to this problem [2][3][4][5]. However associated cost and underdeveloped technology limits the industry from its implementation. Therefore, use of geothermal energy to replace the carbon-based energy sources is gaining momentum [6]. A milestone of the installation of 2 million heat pumps by the European geothermal heat pump market was achieved in 2019 [7]. The geothermal heat usage and electricity production in Europe is expected to grow up to 880-1050 TWh/year and 100-210 TWh/year in 2050 respectively. This contribution is equivalent to 4-7% of European power generation in the year 2050 [8]. As part of the Multidisciplinary and multicontact demonstration of EGS exploration and Exploitation Techniques and potentials [9] project, a numerical hydrothermal model was developed to critically validate the flow behavior of the Soultz-sous-Forêts geothermal power plant from existing operational data. Furthermore, our model was enhanced by including discrete fault structures and validated with operational data to allow for a realistic prediction of the future operational behavior.
Soultz-sous-Forêts is located in the central Upper Rhine Graben, France and has a great potential for geothermal energy exploitation. Soultz-sous-Forêts is the most investigated site in terms of geoscientific studies. The top 1.5 km of the geological succession is made of Soultz-sous-Forêts is located in the central Upper Rhine Graben, France an great potential for geothermal energy exploitation. Soultz-sous-Forêts is the most gated site in terms of geoscientific studies. The top 1.5 km of the geological succe made of thick Quaternary and Tertiary sediments, Mesozoic to Paleozoic sedim rocks above the crystalline basement, which is represented by naturally fractured The Mesozoic to Paleozoic sedimentary rocks can be subdivided into two Buntsandstein and Permian. The Buntsandstein is approximately 350 m thick an prised of fluvial deposits whereas the Permian represents more alluvial continen posit filling the paleo-basin of the variscan orogeny [10]. The basement is comp monzogranite with K-feldspar mega crystals with localized concentration of biotite between 1420 and 4700 m) and a two-mica granite containing muscovite (depth b 4700 and 5000 m) [11,12]. In Table 1, the rock properties for the two sandstone lay granite are listed [13,14]. It must be noted that the data presented in Tables 1 an based on the calibration through the field data and discussed in the unpublished of the MEET project. The sedimentary section has a maximum geothermal gradien to >100 K km −1 making the Soultz-sous-Forêts site ideal for geothermal energy ext [15]. Figure 1 shows that the temperature around the wellbores of Soultz-sous-F higher than that of the surrounding region. Free convection along the major faults is the primary reason causing the increased thermal gradients. For depths great 3700 m, the geothermal gradient becomes 10 K/km.    [13,[20][21][22].
The geothermal project was commenced at Soultz-sous-Forêts in 1984 and the drilling started in 1987 [26]. The earliest plan was to create a fractured granite reservoir in the deep crystalline rock at a depth of 5 km to generate electricity. The industrial electricity production at this site started in June 2016. Presently, the Soultz-sous-Forêts site operates three wells with a maximum depth of up to 5000 m (GPK-2, GPK-3 and GPK-4, see Figure 3). These wells follow the main fault along the NNW-SSE direction. The binary geothermal power plant is working on an Organic Rankine Cycle (ORC) for the heat to electricity conversion. The production well is GPK-2 whereas two wells, GPK-3 and GPK-4, are reinjection wells. The hot fluid produced from GPK-2 is fed into the heat exchanger where the heat is transferred to the isobutane of the ORC cycle and reinjected after being cooled. The fluid production temperature at the Soultz plant is >150 • C and the injection temperature is 70 • C. The production well (GPK-2) and one injection well (GPK-3) indicate fluid leakage in the respective depth intervals at 1431-4170 m measured depth from ground level (MDGL) and 1447-3988 m MDGL, respectively [27]. There is not enough precise data available for the leakage zone and, therefore, it is assumed to be homogeneous over the depth. Both injection wells are cased only at the top, whereas the granitic reservoir section is not completed and in an open-hole condition. the surface to calculate this geothermal gradient. The initial data up to the depth of 5.1 km is measured alongside GPK-2 by Pribnow and Schellschmidt [25] and further modified by Rolin et al. [13].
The geothermal project was commenced at Soultz-sous-Forêts in 1984 and the drilling started in 1987 [26]. The earliest plan was to create a fractured granite reservoir in the deep crystalline rock at a depth of 5 km to generate electricity. The industrial electricity production at this site started in June 2016. Presently, the Soultz-sous-Forêts site operates three wells with a maximum depth of up to 5000 m (GPK-2, GPK-3 and GPK-4, see Figure  3). These wells follow the main fault along the NNW-SSE direction. The binary geothermal power plant is working on an Organic Rankine Cycle (ORC) for the heat to electricity conversion. The production well is GPK-2 whereas two wells, GPK-3 and GPK-4, are reinjection wells. The hot fluid produced from GPK-2 is fed into the heat exchanger where the heat is transferred to the isobutane of the ORC cycle and reinjected after being cooled. The fluid production temperature at the Soultz plant is >150 °C and the injection temperature is 70 °C. The production well (GPK-2) and one injection well (GPK-3) indicate fluid leakage in the respective depth intervals at 1431-4170 m measured depth from ground level (MDGL) and 1447-3988 m MDGL, respectively [27]. There is not enough precise data available for the leakage zone and, therefore, it is assumed to be homogeneous over the depth. Both injection wells are cased only at the top, whereas the granitic reservoir section is not completed and in an open-hole condition.  Table 1 (blue color). Open hole sections of the injection wells are denoted by green colors (GPK-3 and GPK-4) whereas open hole section of the production well is denoted by the dark red color (GPK-2). The leakage zone of the production well is denoted by light red whereas the leakage zone of the GPK-3 is shown by the yellow color.
For the model geometry, only the hydraulically active fractures with high permeability, as proven by thermal anomalies, detected microseismicity during stimulation and operation [23], and which are intersecting multiple wells, were included. The model is thus limited to only five major fractures out of 39 faults or fault zones as shown in Figure 3. The properties of these fractures (fault zones) are listed in Table 2.  Table 1 (blue color). Open hole sections of the injection wells are denoted by green colors (GPK-3 and GPK-4) whereas open hole section of the production well is denoted by the dark red color (GPK-2). The leakage zone of the production well is denoted by light red whereas the leakage zone of the GPK-3 is shown by the yellow color.
For the model geometry, only the hydraulically active fractures with high permeability, as proven by thermal anomalies, detected microseismicity during stimulation and operation [23], and which are intersecting multiple wells, were included. The model is thus limited to only five major fractures out of 39 faults or fault zones as shown in Figure 3. The properties of these fractures (fault zones) are listed in Table 2.
Although the Soultz-sous-Forêts site has been the focus of more than 60 PhD theses and 300 peer-reviewed articles [19], only a few hydrothermal modeling studies have been conducted to understand the hydro-thermal behavior of the reservoir in detail. These studies were coupled with and validated by field operational data specifically with tracer tests to understand the flow path within the fractured granite [14].
The flow circulation between GPK-3 and GPK-2 wells was addressed by Sanjuan et al. [27] through an analytical dispersive transfer model, whereas Blumenthal et al. [28], Gessner et al. [29] and Egert et al. [30] also used dispersive transport models for the Soultz-sous-Forêts site. They investigated the hydraulic connectivity between the injection well (GPK3) and production wells (GPK2 and GPK4) using a multi-well tracer test. Gentier et al. [31] developed the first discrete fracture network (DFN) model while employing a particle tracking method to consider the hydraulically active parts and fracture sets for both wells. More recent modeling studies include Magnenet et al. [32], where a 2D THM model was developed based on a finite element grid (FEM); Aliyu and Chen [33], where finite element method (FEM) was used to model hydro-thermal (HT) processes of Soultz while using different working fluids; and most recently Vallier et al. [14], where a THM model based on FEM was developed at reservoir scale coupled with gravity measurements.
Previous studies showed that a single-fracture approach is not sufficient to represent the hydraulic flow existing at Soultz and 2D models are limited to represent the site in terms of the complex geometry and interconnection of dominating faults. Thus, this study takes its roots from the developed 3D THM model based on FEM while hosting five fractures (FZ1800, FZ2120, FZ4760, FZ4770 and FZ4925; also see Table 2) [13].
From the above literature, it is clear that cold water is injected at 70 • C through both the injection wells. Therefore, injection of cold water below this temperature may enable much higher geothermal energy extraction. However, no numerical studies have been conducted thus far to support this idea. In the presented study, the energy extraction potential from Soultz-Sous-Forêts for 100 years was investigated, allowing the thermal drawdown at the production well to be quantified. The major simplification of this study is neglecting the mechanical behavior. For the short term, as the temperature and pressure development are limited in the wellbore regions, this simplification is relevant and we can use the modeling hydro-thermal simulation result matching with the operational data to better characterize the wellbore effect and reservoir properties. In the ongoing study, we are trying to examine THM behavior of this system for a better prediction for the long term. Another simplification considered here is scaling in the reservoir. Possible scaling effects on the pipelines and heat exchanger devices are beyond the scope of this study. The reservoir size considered for the numerical simulation is large and computational modeling of kinetic controlled reactive fluid flow in such a reservoir requires significantly high computational resources. The possible incompatibility is insignificant because of the reinjection of the same fluid for the entire operation. However, the effect of temperature reduction on the chemical reactions requires experimental work to update the permeability variation.
The manuscript outline is as follows: First, we present a brief geological setting of Soultz-Sous-Forêts, followed by numerical modeling studies for the site. Furthermore, the mathematical and computational technique to model hydro-thermal processes during heat mining from a fractured reservoir is discussed. Next, the wellbore-reservoir coupling is demonstrated and its impact on wellhead temperature is quantified. In the following section, model results and their discussion are followed by final conclusions.

Methodology
In this section, the mathematical modeling is discussed in two stages. In the first part, governing equations for cold water dynamics in the porous media are presented, and in the second part a mathematical model for fluid leakage from the wellbore is discussed.

Reservoir Flow Modeling
A constant heat flux of 0.07 W/m 2 [17]  available, the monthly averaged daily weather fluctuation of Strasbourg, France was used for this study. Strasbourg is approximately 40 km SSE from the Soultz geothermal site. All fractures within the domain are regarded as internal boundaries, implicitly considering the mass and energy exchange between porous media and fractures or fault zones. In the injection well, the diameter of the well is small and can, as a simplification, be represented by a line.
The coupled heat and mass transfer in a fractured rock matrix can be modeled using the mass balance equation integrated with heat transport. The governing equation for heat and mass flow in porous media can be written as [34]: In the above equation, fluid pressure and temperature in the rock matrix are denoted by p and T, respectively. Here, rock porosity is φ m , and storage coefficients for rock and fluid are S 1 and S m . The thermal expansion coefficient of the fluid and rock matrix is denoted by β 1 and β m , respectively. The fluid density and dynamic viscosity are indicated using ρ 1 and µ, whereas the reservoir permeability is denoted by k m .
The fractures are assumed as internal boundaries and the flow along the internal fractures can be denoted by: Here, fluid pressure and temperature in the fracture are indicated by p and T respectively. Additionally, φ f , S f , β f , e h and k f denote the fracture porosity, storage coefficients of the fracture, thermal expansion coefficient of the fracture, hydraulic aperture between the two fracture surfaces, and fracture permeability, respectively. The mass flux exchange between the fracture and matrix are denoted by n.Q m = n.(− ρk m µ∇p ), whereas the gradient operator applicable along the fracture tangential plane is indicated by ∇ T .
The local thermal non-equilibrium (LTNE) approach to model heat exchange between the rock matrix and water is implemented in this study. The conductive heat transfer between rock matrix and pore fluid is the dominant heat exchange mechanism. For the rock matrix, the heat transfer equation can be written as: In the above equation, rock matrix and fluid temperatures are denoted by T m and T l , respectively. Here, rock density, rock-specific heat capacity, rock thermal conductivity and the rock-fluid heat transfer coefficient are denoted by ρ m , C p,m , λ m and q ml , respectively. The heat flux leaving the domain and received by the adjacent fracture can be written as: where T m and T l are the matrix and fluid temperatures in the fracture, respectively; ρ f is the density of the fracture; C p, f is the specific heat capacity of the fracture; λ f is the thermal conductivity of the fracture; and q f l represents the rock fracture-fluid interface heat transfer coefficient, related to the fracture aperture. The last term on the right-hand side of Equation (4) represents the heat flux exchange between the rock matrix and the fracture. The heat convection equation for the pore fluid can be written as: Here C p,l is the heat capacity of the fluid at a constant pressure and λ l is the thermal conductivity of the fluid.
The heat flux coupling relationship of the fluid between the domain and the fracture is satisfied by: where the heat flux n.q l = n.(−φ l λ l ∇T l ) denotes the heat exchange of the fluid between porous media and the fracture. Temperature-dependent fluid thermodynamic properties are implemented into the coupled hydrothermal mass and energy balance equations. The thermophysical properties of water as a function of temperature, including dynamic viscosity (µ), specific heat capacity (C p ), density (ρ) and thermal diffusivity (κ), are listed below [34]: We used the commercial software COMSOL Multiphysics, version 5.6 [34] for numerically solving the coupled mass and energy conservation equations listed above. COM-SOL Multiphysics solves general-purpose partial differential equations using the finite element method.

Wellbore Leakage Modeling
Understanding the fluid flowing temperature along the wellbore can be useful for an accurate estimation of the overall heat production at the production wellhead temperature, and for estimating any possible leakage caused by heat loss along the wellbore. Several reliable analytical techniques are reported in the literature to calculate the flowing temperature distribution along a wellbore [35][36][37].
We integrated our reservoir simulation with a wellbore flow model as developed by Hasan et al. [36]. The model constitutes an analytical approach to estimate wellborefluid temperature distribution for steady state flow. The analytical equations are solved sequentially for each section. Figure 4 shows a simplification of a typical geothermal well with one deviation angle. The well is inclined at an angle α with the horizontal plane. The heat transfer between the wellbore fluid and the rock matrix occurs due to the temperature difference between them. A general energy balance equation for single phase fluid flow can be expressed as:  Here, is the fluid enthalpy, is the gravitational constant, is the depth from the surface, is the flow velocity, is the heat flux per unit and is the mass rate. When assuming no-phase change conditions, enth come: In the above equation, is the fluid temperature and is the press specific heat capacity of fluid and is the Joule-Thomson coefficient. If temperature, the energy balance equation will be: The heat flux per unit wellbore length can be expressed as: Here, is the rock temperature, and is the relaxation parameter In Equations (16) and (17), is the tubing outside radius, is th transfer coefficient, is rock thermal conductivity, is the nondimensi ture, is the measured depth of the wellbore, is the geothermal grad Here, H is the fluid enthalpy, g is the gravitational constant, z is the variable well depth from the surface, ν is the flow velocity, Q is the heat flux per unit of well length and w is the mass rate. When assuming no-phase change conditions, enthalpy will become: In the above equation, T is the fluid temperature and p is the pressure, c p is the specific heat capacity of fluid and C J is the Joule-Thomson coefficient. If T f is the fluid temperature, the energy balance equation will be: The heat flux per unit wellbore length can be expressed as: Here, T ei is the rock temperature, and L R is the relaxation parameter defined as: In Equations (16) and (17), r to is the tubing outside radius, U to is the overall heat transfer coefficient, k e is rock thermal conductivity, T D is the nondimensional temperature, L is the measured depth of the wellbore, g G is the geothermal gradient and Φ is the lumped parameter, which lumps the kinetic energy term and the Joule-Thomson coefficient term.
If V is the fluid specific volume and S is fluid entropy then from Maxwell identities, we can write: For liquids where ρ is the liquid density, volume expansivity (β) can be calculated as: Therefore, the final output temperature from the wellhead will be: In this text, we considered three wells: GPK-3 and GPK-4 as two injection wells and GPK-2 as a production well. In GPK-3, the wellbore leakage was assumed between 1282 and 4852 m depth measured from the surface. In the case of GPK-2, the wellbore leakage was modeled between 1264 m to 4244 m depth measured from the surface. The fluid is single phase water flow and the model parameters are constant specific heat capacity of water as 4200 J·kg −1 K −1 , L R = 0.00001 m −1 , and Φ = 0.00345 Km −1 , respectively. Here, L R and Φ accounts for the casing properties, cement properties and their thicknesses.
The coupling between the reservoir and the wellbore model is achieved through a sequential approach. First, the temperature drop due to heat exchange between the injection wellbore and the rock matrix is calculated through the analytical model. From this, the final wellbore bottom temperature is obtained, which is used as an input for the first iteration of the numerical reservoir model (heat exchange between the rock and the fluid). In the next stage, the wellbore heat exchange effect is implemented through the updated values for the reservoir temperature measured at the production wellbore bottom. The wellbore heat exchange effect is defined analytically and the temperature alongside the wellbore is obtained. Wellbore radius is very small compared to the reservoir size, and is considered as a line with the calculated temperature profile through the analytical model inside the reservoir simulation. The total number of elements in the geometry is 142,051, whereas boundary elements number 8305 and edge elements number 666.

Results and Discussions
In this section, first we present the benchmark for our numerical model. Then, the hydrothermal numerical modeling results are compared with the operational data measured at Soultz-sous-Forêts for three years of operation. Furthermore, new injection scenarios are proposed that can be adopted with the existing industrial setup to enhance the energy extraction capability. Finally, we perform sensitivity analysis on ten governing parameters and estimate their impact on the production temperature.

Benchmarking
For benchmarking the numerical model, we used the approach adopted by Cheng et al. [38] and Bongole et al. [39] by using a simplified 1D heat transfer problem for a single fracture system. This approach is used in the previous studies to benchmark the models. The analytical equation for heat transfer considers that the geometry is infinitely extended in both directions (see Figure 5), there is no flow boundary conditions for heat exchange, steady state fluid flow occurs only through the fracture and the rock permeability is zero, and the thermophysical properties of water are constant throughout the simulation. The temperature distribution for the fluid is identical to that of the rock matrix due to the local thermal equilibrium assumption. The analytical solution for the fluid temperature distribution is [38,39]: thermal equilibrium assumption. The analytical solution for the fluid temp bution is [38,39]:  We observe a good agreement between the numerical and analytica demonstrated in Figure 6. The operational data for three years was made available for Soultz-s by the site operators and is used here to calibrate the coupled unsteady model. Figure 7 shows the injection and production rates at the wellhead from June 2016 to September 2019. The fluid injection temperature is 70 ° injection wells. In the above equation, T f D is the nondimensional fracture temperature and u f is the fluid velocity.
We observe a good agreement between the numerical and analytical solutions, as demonstrated in Figure 6. thermal equilibrium assumption. The analytical solution for the fluid temperature distribution is [38,39]:

( − ) ]
In the above equation, is the nondimensional fracture temperature and is the fluid velocity. We observe a good agreement between the numerical and analytical solutions, as demonstrated in Figure 6. The operational data for three years was made available for Soultz-sous-Forêts site by the site operators and is used here to calibrate the coupled unsteady hydro-thermal model. Figure 7 shows the injection and production rates at the wellhead for 1163 days from June 2016 to September 2019. The fluid injection temperature is 70 °C for both the injection wells. The operational data for three years was made available for Soultz-sous-Forêts site by the site operators and is used here to calibrate the coupled unsteady hydro-thermal model. Figure 7 shows the injection and production rates at the wellhead for 1163 days from June 2016 to September 2019. The fluid injection temperature is 70 • C for both the injection wells.

Validation with Operational Data
In Figure 8, the numerical model data is validated with operational data for the tim period as described above. Unfortunately, it is not possible to publish the exact values o operational data due to concerns of our industrial partners and we can only show th amount of change. These are the actual temperature values shown by different colors fo operational and simulations (not the differences). The measured temperature data is th operational data for 1163 days at the production wellhead. The temperature at the pro duction well based on the hydro-thermal model is significantly different compared to th operational data. For most of the operational period, the predicted production well tem perature is 15 °C higher than the measured temperature. Only operation onset and term nation stages display smaller deviation in predicted temperature than observed tempera ture. Because the wellhead temperature measuring device may be affected by the loc ambient temperature and the monthly average temperature near the geothermal site almost the same for the corresponding months in each operational year, a correction facto

Validation with Operational Data
In Figure 8, the numerical model data is validated with operational data for the time period as described above. Unfortunately, it is not possible to publish the exact values of operational data due to concerns of our industrial partners and we can only show the amount of change. These are the actual temperature values shown by different colors for operational and simulations (not the differences). The measured temperature data is the operational data for 1163 days at the production wellhead. The temperature at the production well based on the hydro-thermal model is significantly different compared to the operational data. For most of the operational period, the predicted production well temperature is 15 • C higher than the measured temperature. Only operation onset and termination stages display smaller deviation in predicted temperature than observed temperature. Because the wellhead temperature measuring device may be affected by the local ambient temperature and the monthly average temperature near the geothermal site is almost the same for the corresponding months in each operational year, a correction factor to account for the weather impact on the measuring device based on the numerical model is introduced. Two scenarios of seasonal impact on the production fluid temperature are considered: (a) 20% impact of ambient temperature (T effective =T simulation + 0.2 × ambient temperature) and (b) 50% impact of ambient temperature (T effective =T simulation + 0.5 × ambient temperature).
iences 2021, 11,464 is introduced. Two scenarios of seasonal impact on the productio considered: (a) 20% impact of ambient temperature (Teffective=Tsimu perature) and (b) 50% impact of ambient temperature (Teffective = temperature).  Figure 8 shows the comparison of the operational data wit wellbore model and the weather-influenced production fluid tem wellbore-reservoir model has the highest overestimation of p However, when daily weather fluctuations in the integrated wellb considered, the prediction matches very well for most of the opera 8. The temperature differences are more relevant to understand i tual value rather than the original temperature. The difference b numerical data while considering 50% of the ambient temperatur perature of the coupled wellbore-reservoir model has the best m els. However, the model deviates by more than 15 °C from the op periods of 1.8 and 2.4 years. Because no other reasons for these with the operational data set, different measurement procedure at the wellhead are assumed as reasons for these deviations.

Long-Term Operational Behavior
In the next study, the model was extended to a simulation operation to predict the wellhead temperature development at the section, different initial temperatures at the bottom hole section measured data were used. The main objective of this study was  Figure 8 shows the comparison of the operational data with the coupled reservoirwellbore model and the weather-influenced production fluid temperature. The integrated wellbore-reservoir model has the highest overestimation of production temperature. However, when daily weather fluctuations in the integrated wellbore-reservoir model are considered, the prediction matches very well for most of the operation, as shown in Figure 8. The temperature differences are more relevant to understand its deviation from the actual value rather than the original temperature. The difference between operational and numerical data while considering 50% of the ambient temperature on the production temperature of the coupled wellbore-reservoir model has the best matching among all models. However, the model deviates by more than 15 • C from the operation data during the periods of 1.8 and 2.4 years. Because no other reasons for these deviations are provided with the operational data set, different measurement procedures or false measurements at the wellhead are assumed as reasons for these deviations.

Long-Term Operational Behavior
In the next study, the model was extended to a simulation period of 100 years of operation to predict the wellhead temperature development at the production well. In this section, different initial temperatures at the bottom hole section than the operationally measured data were used. The main objective of this study was to estimate the temperature at the production well (GPK-2) for different injection temperatures for long-term operational periods. In both scenarios, the injection rates for the first 1163 days are the same as in the provided operational data set. The recently designed heat exchanger at Soultz-sous-Forêts is capable of cooling the water from 70 to 40 • C using a cooling loop at 15 • C and 40 m 3 /h [40]. Therefore, two scenarios were considered, A and B, for different injection temperatures. For the remaining operational period, scenario A considers four different fluid injection temperatures at the injection wellhead (70, 60, 50 and 40 • C).
The fluid injection rates are 13.3 and 11 L/s for GPK-3 and GPK-4, respectively, and the production fluid rate of GPK-2 is 24.3 L/s for the remaining operational period. In Scenario B, the injection rates after 1163 days are 19.6 and 9.7 L/s for GPK-3 and GPK-4, respectively, and the production rate of GPK-2 is 29.3 L/s; the same four injection wellhead temperatures as for scenario A were considered: 70, 60, 50 and 40 • C. These values of the injection and production rates are the operational requirements requested by our industrial partner. Figure 9a-d shows the temperature along the wellbore for scenarios A and B, respectively, for both injection wells. The wellbore GPK-3 has an open hole section that causes a linear temperature drop along the wellbore instead of a nonlinear temperature drop, as shown in Figure 9a,b. It is interesting to note that instead of having different injection-production rates in all three wells, the fluid production temperature at the GPK-2 wellhead is almost similar for both of the scenarios A and B, as shown in Figure 9e,f. The small increase in temperature at the production wellhead is due to the sudden drop in the production wellhead pressure. The contribution in the fluid flow is due to the first pressure shock of the injection that comes from the faulted zones which are located at the bottom of the system with a higher temperature. As time proceeds, the contribution from the matrix and the leakage zone increases and reduces the temperature a few days after the beginning of the injection. To calculate the initial temperature at the wellhead, it is assumed that there is a steady state flow from the combination of the matrix and the fault zones. This initial temperature is slightly lower than that of the unsteady condition at the early time period. The fluid with the lower viscosity shows a delay in the development of the pressure shock resulting from the cold fluid injection. Therefore, the contribution from the matrix and the leakage zone for the fluid with the temperature 40 • C happens later and the main fluid flow from the faulted zone in the bottom of the system lasts for a longer time. Moreover, the temperature increase in scenario B is higher compared to that of scenario A due to the fact that scenario B has a higher production rate than scenario A which reduces the time for exchanging heat in the wellbore. Figure 10 shows the comparison of temperature distribution in the fractures and along the wellbore for scenarios A and B. The higher production rate results in slightly faster thermal drawdown at the production well bottom for scenario B than scenario A. No thermal breakthrough was observed at the production well bottom even after 100 years of operation, as shown in Figure 10e,f.

Uncertainties
There are several uncertainties in this model. We considered the wellbore as a line source for the heat flow. The faulted zone is formulated using a fracture. Both of these assumptions are reliable because the size of the wellbore and the fault zone is negligible in comparison to the overall size of the reservoir. Data validation for the short operational period for production confirms this behavior. The matrix zone is considered as homogeneous and isotropic. As the permeability of matrix is lower than faulted zone, its contribution to the heat and mass flux is small. Therefore, this assumption holds true. As we do not know the exact point of the leakage zone alongside the casing area, we considered a homogeneous leakage and tried to compensate for the possible errors by performing a trial-and-error method to find an appropriate lumped parameter that defines the wellbore heat exchange effect. Furthermore, due to the unavailability of the geomechanical and geochemical data, we mainly focused on the hydrothermal behavior of the geothermal system. Short-term validation of this TH model gives an insight regarding the accurate system characterization, including the permeability and porosity distribution, fault placement and its contribution to the overall flow, the wellbore effect on the overall heat exchange, and fluid and rock properties. Therefore, it builds a basis for future THM or THMC (thermo-hydro-mechanical-chemical) simulations. Our expectation regarding the THM behavior is that permeability around the injection would increase (resulting from the localized thermoelastic stress reorientation and increased pore pressure). Therefore, the enhanced permeability will be favorable in the energy extraction. Based on this un-derstanding, TH provides a preliminary basis for the thermoelastic and poroelastic stress development in this geothermal site.
This initial temperature is slightly lower than that of the unsteady condition at the ea time period. The fluid with the lower viscosity shows a delay in the development of pressure shock resulting from the cold fluid injection. Therefore, the contribution from matrix and the leakage zone for the fluid with the temperature 40 °C happens later a the main fluid flow from the faulted zone in the bottom of the system lasts for a long time. Moreover, the temperature increase in scenario B is higher compared to that of s nario A due to the fact that scenario B has a higher production rate than scenario A wh reduces the time for exchanging heat in the wellbore.  Geosciences 2021, 11,464 14 of 19 Figure 10 shows the comparison of temperature distribution in the fractures and along the wellbore for scenarios A and B. The higher production rate results in slightly faster thermal drawdown at the production well bottom for scenario B than scenario A. No thermal breakthrough was observed at the production well bottom even after 100 years of operation, as shown in Figure 10e,f.

Uncertainties
There are several uncertainties in this model. We considered the wellbore as a line source for the heat flow. The faulted zone is formulated using a fracture. Both of these assumptions are reliable because the size of the wellbore and the fault zone is negligible in comparison to the overall size of the reservoir. Data validation for the short operational period for production confirms this behavior. The matrix zone is considered as homogeneous and isotropic. As the permeability of matrix is lower than faulted zone, its contribution to the heat and mass flux is small. Therefore, this assumption holds true. As we do not know the exact point of the leakage zone alongside the casing area, we considered a homogeneous leakage and tried to compensate for the possible errors by performing a trial-and-error method to find an appropriate lumped parameter that defines the wellbore heat exchange effect. Furthermore, due to the unavailability of the geomechanical and geochemical data, we mainly focused on the hydrothermal behavior of the geothermal system. Short-term validation of this TH model gives an insight regarding the accurate system characterization, including the permeability and porosity distribution, fault placement and its contribution to the overall flow, the wellbore effect on the overall heat exchange, and fluid and rock properties. Therefore, it builds a basis for future THM or THMC (thermo-hydro-mechanical-chemical) simulations. Our expectation regarding the

Sensitivity Analysis of Hydrothermal Uncertainties
To examine the effect of the uncertainty of the involved parameters in hydro-thermal simulations for the Soultz-sous-Forêts geothermal reservoir, a detailed sensitivity analysis was performed, and its results are shown in Figure 11. The base case was selected as scenario A with an injection temperature of 40 • C. The range of parameters are mentioned in Table 3. Results show that the hydraulic conductivity of the matrix and fault zone and the wellbore leakage have a considerable effect on the production temperature. This finding is well aligned with the sensitivity analysis of the THM process in the fracture reservoir system [41,42]. However, these values are one order of magnitude lower than the wellbore heat exchange effect, as shown in Figure 9e. The other three parameters-fault thickness, matrix thermal conductivity and the matrix specific heat capacity-have approximately 1 • C variation over 100 years of operation. The porosity of the matrix and the fault zone, in addition to the fault zone thermal conductivity, have no impact on the temperature variation. Interestingly, heat flux has no effect on the production temperature at the surface due to the conductive heat flow in the reservoir and because the fault zones are farther away from the bottom boundary considered for the simulation. Important parameters in this sensitivity analysis show a monotonic effect on the production temperature behavior and they cannot explain the sinusoidal temperature fluctuation of more than 10 • C in each year. To check our assumption regarding the weather fluctuation impact on the production temperature, we tried to estimate the wellbore heat exchange effect. We found that the wellbore heat exchange is mainly flow-rate dependent parameter and the flow rates for the production data are constant from 41 to 250 days, whereas we can see a fluctuation in the recorded temperature. This indicates that the cyclic variability in the production temperature cannot be supported by the wellbore heat exchange argument. Therefore, the most suitable reason of the periodic production temperature variability is weather fluctuation on the measuring device. Figure 11. Sensitivity analysis for 10 parameters affecting the hydro-thermal processes at Soultz-sous-Forêts for (a) matrix hydraulic conductivity, (b) heat flux from the bottom boundary, (c) matrix specific heat capacity, (d) hydraulic conductivity of faults (here , is the fault zone hydraulic conductivity as given in Table 2), (e) porosity of fault zone, (f) leakage contribution to the total fluid flow, (g) matrix porosity, (h) matrix thermal conductivity, (i) fault thickness (here is the fault thickness as given in Table 2), and (j) thermal conductivity of the fault zone.

Conclusions
As part of the MEET project, a coupled reservoir and wellbore model for hydraulic and thermal processes involved during the geothermal energy extraction operation at Soultz Sous Forêts was developed. Operational data from a period of 1163 days of operation was used to validate the numerical model. The validated hydro-thermal numerical model precisely simulates the geothermal energy extraction operation for 3 years. Furthermore, two operational scenarios for 100 years with four different injection wellhead temperatures-70, 60, 50 and 40 °C-were analyzed. It can be observed that even after 100 years of operation, the thermal breakthrough at the production well is only in the range of 10 to 20 °C. After 100 years of cold fluid injection and hot fluid production, the observed temperature drop at the production wellhead is less than 20 °C . Therefore, our numerical model predicts that 100 years of geothermal energy extraction operation at Soultz-sous- Figure 11. Sensitivity analysis for 10 parameters affecting the hydro-thermal processes at Soultz-sous-Forêts for (a) matrix hydraulic conductivity, (b) heat flux from the bottom boundary, (c) matrix specific heat capacity, (d) hydraulic conductivity of faults (here K f ,0 is the fault zone hydraulic conductivity as given in Table 2), (e) porosity of fault zone, (f) leakage contribution to the total fluid flow, (g) matrix porosity, (h) matrix thermal conductivity, (i) fault thickness (here F 0 is the fault thickness as given in Table 2), and (j) thermal conductivity of the fault zone. Hydraulic conductivity of fault zone (see Table 2)   Table 2) F 0 m 0.5F 0 m 2F 0 m Thermal conductivity of the fault zone 2.5 W/m/K 2 W/m/K 3 W/m/K

Conclusions
As part of the MEET project, a coupled reservoir and wellbore model for hydraulic and thermal processes involved during the geothermal energy extraction operation at Soultz Sous Forêts was developed. Operational data from a period of 1163 days of operation was used to validate the numerical model. The validated hydro-thermal numerical model precisely simulates the geothermal energy extraction operation for 3 years. Furthermore, two operational scenarios for 100 years with four different injection wellhead temperatures-70, 60, 50 and 40 • C-were analyzed. It can be observed that even after 100 years of operation, the thermal breakthrough at the production well is only in the range of 10 to 20 • C. After 100 years of cold fluid injection and hot fluid production, the observed temperature drop at the production wellhead is less than 20 • C. Therefore, our numerical model predicts that 100 years of geothermal energy extraction operation at Soultz-sous-Forêts is feasible and will have a sufficiently high production temperature throughout the operation duration.