# Large-scale reservoir modeling of the Vendenheim geothermal site (France)

Javier Abreu-Torres, Gergő Hutka, Guido Blöcher, Mauro Cacace, Vincent Magnenet, Jean Schmittbuhl

► To cite this version:

Javier Abreu-Torres, Gergő Hutka, Guido Blöcher, Mauro Cacace, Vincent Magnenet, et al. Large-scale reservoir modeling of the Vendenheim geothermal site (France). Advances in Geosciences, 2025, 65, pp.117 - 125. 10.5194/adgeo-65-117-2025. hal-04998701

HAL Id: hal-04998701

https://hal.science/hal-04998701v1

Submitted on 20 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Large-scale reservoir modeling of the Vendenheim geothermal site (France)

Javier Abreu-Torres<fnref n="1" />, Gergő Hutka<fnref n="2" />,<fnref n="3" />, Guido Blöcher<fnref n="2" />,<fnref n="3" />, Mauro Cacace<fnref n="2" />, Vincent Magnenet<fnref n="1" />, and Jean Schmittbuhl<fnref n="1" />

**Correspondence:** Javier Abreu-Torres (abreutorres@unistra.fr)

Received: 30 June 2024 – Revised: 4 November 2024 – Accepted: 8 November 2024 – Published: 6 January 2025

**Abstract.** During the Vendenheim deep geothermal project (Strasbourg Eurometropole, France), large induced seismic events led to the arrest of the project. Two important features of the induced seismicity were unexpected: the large distance to the wells of a cluster of seismic events (4–5 km) and the occurrence of the largest event Mlv3. 9 at the bottom of the wells, six months after shut-in. To better understand the natural hydro-thermal conditions at the Vendenheim site, we develop within the framework of the DT-GEO project (Horizon Europe) a simplified large-scale model (8 km × 8 km × 6 km) of the area. We aim at performing in-silico experimentation to reproduce the natural hydrothermal state of the geothermal reservoir related to the heat flow (conductive/convective) within the model. We first test our methods by solving 2D and 3D benchmarks related to the convective regime in saturated porous media. Our model is based on the MOOSE/-GOLEM framework (finite element approach) and integrate the public regional geological model GEORG that includes major lithologies. We present the coarse-grained simulations of the natural fluid circulation.

## 1 Introduction

The Vendenheim geothermal site is situated in the Upper Rhine Graben region (URG), which is located in the Rhine valley between the Vosges and the Black Forest mountains, in borders of France and Germany, and in the north of Switzerland. The URG is one of the promising regions for geothermal exploration in the European context (Frey et al., 2022).

However, it has been studied and verified that in the context of geothermal systems, induced seismicity could be generated, causing many projects to be abandoned, (Majer et al., 2007; Zang et al., 2014; Kwiatek et al., 2019; Gaucher et al., 2015; Schmittbuhl et al., 2021). This was the case of the Geoven project that was developed in the Vendenheim area, in the north of Strasbourg France. This project was closed due to a series of earthquakes between March 2018 and June 2021, in which the largest event had a magnitude of $M_{LV}$ 3.6 (Schmittbuhl et al., 2021).

For this reason, as part of the Digital Twin for GEOphysical Extremes (DT-GEO) project (https://dtgeo.eu, last access: 4 October 2024), we develop a 3D Thermal Hydrological coupled (TH) modelling to compute the hydrothermal flow with a temperature dependent density and viscosity rheology, in a porous medium, using the MOOSE/GOLEM application (Jacquey and Cacace, 2017; Giudicelli et al., 2024), in order to verify if a hydrothermal convective system is compatible with known observations at the Vendenheim site. Furthermore, we aim to get a better idea of the initial conditions of a model for induced seismicity, as this hydrothermal convection could introduce characteristic time and space scales.

We use a simplified structural model based on the GEORG platform (https://maps.geopotenziale.eu, last access: 8 October 2024), within the Vendenheim area (Strasbourg, France) in order to set the distribution of the principal geological layers. In the shallow crust fluid circulation in geological systems, such as: tectonically active systems, ash-flow calderas, sedimentary basins and fault regions, have been shown to exist. These systems may have a rock permeability large

<footnote n="1">¹Ecole et Observatoire des Sciences de la Terre, Strasbourg University/CRNS, 5 Rue René Descartes, 67000 Strasbourg, France</footnote>

<footnote n="2">²German Research Centre for Geosciences, Helmholtz-Zentrum Potsdam, Deutsches GeoForschungsZentrum GFZ, 14473 Potsdam, Germany</footnote>

<footnote n="3">³TU Berlin, Department of Engineering Geology, 10587 Berlin, Germany</footnote>
enough to allow thermally driven fluid flow (Davaille, 1999; Simms and Garven, 2004; Lopez et al., 2016; Patterson et al., 2018; Duwiquet et al., 2019; Guillou-Frottier et al., 2020). Thermally driven fluid circulation can occur in a variety of systems in Earth's interior (Donaldson, 1970; Schubert et al., 2001; Bataillé et al., 2006; Angiboust et al., 2012).

In the context of geothermal fields as Vendenheim, the geothermal potential of an area, which is based on environmental, technical, economic and acceptance principle according to Bayer et al. (2019), must be estimated by evaluating the temperature field and the mode of heat transfer (Ranalli and Rybach, 2005). In a saturated porous medium heated from below (i.e. higher temperature at depth than at surface), heat transfer is better understood by the Rayleigh number, that indicates the relative importance of the viscous stabilising mechanisms to the buoyancy destabilising mechanisms. The Rayleigh number in porous media is defined according to Schubert and Straus (1979), Nield (1994), Desaive (2003) as:

$$Ra = \frac{k \rho_f c_f g \beta (T_b - T_t) H}{\lambda_b \mu}, \quad (1)$$

where k is the intrinsic permeability, $\rho_f$ is the density of the fluid, $c_f$ is the fluid heat capacity, g is the gravity acceleration, $\beta$ is the fluid thermal expansion, $T_b$ and $T_t$ are the bottom and top temperature respectively, H is the height of the model, $\lambda_b$ is the bulk heat conductivity, and $\mu$ is the fluid viscosity. The critical Rayleigh number ($Ra_c$) is defined as the value of the Rayleigh number at which instability begins and free convection develops: When $Ra > Ra_c$, the system is unstable (heat is transferred by convection), and when $Ra < Ra_c$, the system is stable (heat is transferred by conduction. In porous media we find in the literature that $Ra_c = 4\pi^2$ (Lapwood, 1948; Murphy, 1979; Desaive, 2003). Furthermore, thermal convection is controlled by two main processes: an increased pressure gradient (“forced convection”), and buoyancy effects due to temperature-dependent fluid density.

In this context, we found Pola et al. (2020), in which a representative model from NE Italy with low-temperature gradient, was developed, in which thermal convection was observed. Moreover, in Liu et al. (2021), a 3D model with a fracture the type of heat flow was investigated using the Nusselt number. Furthermore in Detournay et al. (2022), they developed a 3D fracture model, where 3 heat transfer mechanisms were considered in a geothermal context. In contrast to these works, here we discuss the Rayleigh number to better understand the heat transfer in the Vendenheim region.

In the following sections we present our methods for TH modeling (Sect. 2). After that, in Sect. 3 we present 2D and 3D benchmarks to test our methods in the context of convective flow. In Sect. 4 we present the simplified model based on the Vendenheim site, with the results of the TH modeling, and discuss the possibility of large-scale natural hydrothermal circulation in the geothermal reservoir. The mechanical consequences of possible thermal anomalies (fault instability) as well as numerical stability (mesh sensitivity) require further analysis and are beyond the scope of this article.

## 2 Methods

In the following, the governing equations are present, in the context of the Thermal Hydrological coupling (TH). This system is solved by Golem (Jacquey and Cacace, 2017; Cacace and Jacquey, 2017), a Moose-based application (Guedicelli et al., 2024) for modelling coupled-Thermo-Hydro-Mechanical processes in faulted geothermal reservoirs.

First, the equilibrium of fluid mass equation in a non-deformable porous matrix is presented:

$$\frac{n}{K_f} \frac{\partial p_f}{\partial t} + \nabla \cdot \mathbf{q}_D = 0, \quad (2)$$

where n is the porosity, $K_f$ represents the fluid bulk modulus, t is the time, and $\mathbf{q}_D$ is the Darcy velocity, which is described according to the Darcy's law as:

$$\mathbf{q}_D = -\frac{k}{\mu_f} \cdot (\nabla p_f - \rho_f \mathbf{g}), \quad (3)$$

where $k$, $\mu_f$, $\rho_f$, and $g$ represent the permeability tensor, the viscosity and fluid density, and the gravity acceleration respectively.

The temperature distribution is governed by the energy balance equation:

$$(\rho c)_b \frac{\partial T}{\partial t} + \nabla \cdot \mathbf{q}_T = 0, \quad (4)$$

where $(\rho c)_b$ is the product of the bulk density and heat capacity, and $\mathbf{q}_T$ is the heat flux, defined as:

$$\mathbf{q}_T = \rho_f c_f \mathbf{q}_D T - \lambda_b \nabla T, \quad (5)$$

where $c_f$ is the specific heat capacity of the fluid and $\lambda_b = n\lambda_f + (1-n)\lambda_s$ is the bulk thermal conductivity, with $\lambda_f$ and $\lambda_s$ being the fluid and solid thermal conductivity. The first and second terms correspond to the convective and conductive flux respectively. For further information, we refer to Magnenet et al. (2014), Blöcher et al. (2018), Vallier et al. (2019).

## 3 2D and 3D thermal convection benchmarks in faulted systems

To validate our methods, we reproduce the benchmarks proposed in Magri et al. (2017) and Guillou-Frottier et al. (2020), where the Elder problem is solved in the context of 2D and 3D hydrothermal convection in fault zones. We restrict our study to the hydrothermal problem, no mechanics are considered.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Temperature distribution (a), and horizontal temperature profile (b) for $k = 1 \times 10^{-15} m^2$ at a depth of 0.75 km. The profile passes through the white line shown in panel (a). The steady state is observed in this example using a tolerance of $1 \times 10^{-4}$. (See Fig. 1 in Guillou-Frottier et al. (2020). The model is 4.5 km long and 1.5 km in depth.</figcaption></figure>

<table><caption>**Table 1.** Boundary conditions used in Guillou-Frottier et al. (2020) for the 2D results.</caption><thead><tr><th>Variable</th><th>Top</th><th>Bottom</th><th>West</th><th>East</th></tr></thead><tbody><tr><td>Pore pressure</td><td>$p_f = 1 \times 10^5$ Pa</td><td>$\nabla p_f = 0$</td><td>$\nabla p_f = 0$</td><td>$\nabla p_f = 0$</td></tr><tr><td>Temperature</td><td>$T = 0$ °C</td><td>$T = 400$ °C</td><td>$\nabla T = 0$</td><td>$\nabla T = 0$</td></tr></tbody></table>

<table><caption>**Table 2.** Values of the physical parameters used to evaluate the convective solution. Porosity values are on a scale from 0 to 1. Permeabilities are isotropic.</caption><thead><tr><th>Solid Properties</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Porosity <i>n</i></td><td>[1]</td><td>0.2</td></tr><tr><td>Solid density <i>ρ</i><sub>s</sub></td><td>kg m<sup>-3</sup></td><td>2700<sub>1</sub></td></tr><tr><td>Specific solid heat capacity <i>c</i><sub>s</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>800.0</td></tr><tr><td>Solid conductivity <i>λ</i><sub>s</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>3.0</td></tr><tr><td>Isotropic permeability <i>k</i></td><td>m<sup>2</sup></td><td>1 × 10<sup>-15</sup></td></tr><tr><td colspan="3"><strong>Fluid properties</strong></td></tr><tr><td>Top fluid density <i>ρ</i><sub>f0</sub></td><td>kg m<sup>-3</sup></td><td>1036.5</td></tr><tr><td>Top fluid viscosity <i>μ</i><sub>f0</sub></td><td>Pa s</td><td>0.001</td></tr><tr><td>Fluid conductivity <i>λ</i><sub>f</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>0.65</td></tr><tr><td>Specific fluid heat capacity <i>c</i><sub>f</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>4180</td></tr></tbody></table>

## 3.1 2D benchmark

The model dimensions are 1.5 km-thick (z direction) and 4.5 km-wide (x direction), with square elements ($Δz = Δx = 20$ m), having 16 875 elements in total. Fixed temperatures of 0 and 400 °C are imposed in the top and bottom boundary respectively, no temperature flow is considered at the sides. For the pore pressure, no flow is imposed in all boundaries except the top boundary. Here a fixed pressure of $1 \times 10^5$ Pa is imposed. A summary of these conditions is presented in Table 1.

The material properties are reported in Table 2.

Fluid density $ρ_f$ and fluid viscosity $μ_f$ obey the following laws:

$$ \begin{aligned} \rho_f(T) &= \rho_{f0} - 0.14167 \times T - 0.0022381 \times T^2 \\ \mu_f(T) &= \mu_{f0} \times 10^{\frac{247.8}{T+133}}, \end{aligned} \quad (6) $$

with $T$ in °C, $ρ_{f0} = 1036.6$ kg m$^{-3}$, $μ_{f0} = 2.414 \times 10^{-5}$ Pa s, are the fluid density and viscosity at the top of the model.

We implement these laws using a Golem function that allows us to introduce the variation of density and viscosity with temperature using tables. These tables are elaborated externally by applying the function in the desired temperature range.

We impose initial conditions of $T(t=0) = (T_b - T_t)/1500 \cdot z$ (with $T_b$ and $T_t$ being the temperature at the bottom and top of the model), and $p_f(t=0) = ρgz$. We evaluate the convergence of the system towards the steady state by using the following convergence criteria:

$$ \frac{u^{n+1} - u^n}{\Delta t} < \tau u^{n+1}, \quad (7) $$

where $u$ represents the solution of the system, $Δt$ is the time step, the superscript $n+1$ and $n$ refer to the current and old solutions respectively, and $τ$ is the tolerance ($τ = 1 \times 10^{-14}$ in this case). This equation represents definition of the L2-norm on relative change in the solution between two time steps.

In Fig. 1, we reproduce the 2D results of Guillou-Frottier et al. (2020) (Fig. 1 in that paper) for the permeability $k = 1 \times 10^{-15} m^2$, in 91.958 kyr at a depth of 0.75 km. The horizontal profile of temperature in the middle of the model
<figure><img src="image_4.png" /><figcaption>Figure 2. Variation of the maximum value of the vertical component of Darcy's velocity with respect to the permeability from $5 \times 10^{-16}$ to $1 \times 10^{-13}$ m$^2$. Red marks represent the permeability values where convection develops, while for blue marks values no convection is observed. Both axes are in logarithmic scale. In the convective regime (red dots), the behavior is fitted by a linear relationship: $v_D = m \cdot k + b$, where m is $1.61 \times 10^7$ and b is $-8.33 \times 10^{-9}$.</figcaption></figure>

(showing by the white line) compared to profile obtained in Guillou-Frottier et al. (2020). The average temperature $T_{avg}$ (232 °C in this work and 236 °C in theirs), and $\Delta T$, defined as the difference between the maximal and average temperature within the profile, are also very similar (83 °C in this work and 84 °C in theirs). The steady state is reached in our work with a tolerance of $1 \times 10^{-14}$.

### 3.1.1 Darcy velocity as a function of the permeability

Keeping the same model, we studied the behavior of the maximum Darcy velocity with respect to different isotropic permeabilities in a range from $5 \times 10^{-16}$ to $1 \times 10^{-13}$ m$^2$. In Fig. 2 we show the variation of the maximal value of the vertical component of the Darcy velocity with respect to the permeability. The red marks represent the permeability values where convection is observed, while blue marks show the permeability values where convection is not observed. For permeability values where convection develops, we select the maximum Darcy velocity value over time. For permeabilities of $5 \times 10^{-16}$ and $6.25 \times 10^{-16}$ m$^2$, convection is not observed, therefore Darcy velocity is very low. We observe that the Darcy velocity increases in the convective regime when the permeability increases. Moreover, for permeability values for which convection is observed, the relationship between permeability and Darcy velocity appears to be linear. This is demonstrated by the dashed line in Fig. 2, which illustrates this relationship and was obtained through linear regression. This can also be explained by the Eq. (3), which states that the Darcy velocity is directly proportional to the permeability if the pressure difference is constant.

<figure><img src="image_5.png" /><figcaption>Figure 3. Faulted model used for the 3D benchmark according to Magri et al. (2017). The model has 5.5 km in the x, y, z directions, and contains a vertical fault with 40 m wide, and 5.5 km long, located at the center of the model, parallel to the y axis. The hexahedrons are 40.14 m on each side in the enclosing rock, and 10.03 m on each side within the fault.</figcaption></figure>

## 3.2 3D benchmark results

The geometry of the model for the 3D benchmark, proposed in Magri et al. (2017), is composed by a 40 m wide, 5.5 km long fault, surrounded by a rock matrix ($5.5 \times 5.5 \times 5.5$) km. The vertical fault is located in the center of the model, parallel to the y axis. In Fig. 3 we show the geometry of the model. Cubic meshes are used, the meshing within the fault is finer. The cubes are 40.14 m on each side in the enclosing rock, and in the fault they are 10.03 m on each side.

In consideration of the boundary conditions, the system is heated with a constant temperature of $T_b = 170$ °C at the bottom of the model, while a constant temperature of $T_c = 20$ °C is imposed at the top of the model. It is assumed that no thermal flow occurs on the remaining sides of the model. A pressure of $1 \times 10^5$ Pa is imposed at the top of the model, with no flow on the other sides. These boundary conditions are expressed in Table 3. We report the properties of the enclosing rock, fault, and fluid in Table 4.

The fluid density $\rho_f$ and viscosity $\mu_f$ follow the next laws:

$$ \begin{aligned} \rho_f(T) &= \rho_{f0} \cdot (1 - \beta(T - T_t)) \\ \mu_f(T) &= \mu_{f0} \cdot \exp\left(-\frac{T - T_t}{T_v}\right), \end{aligned} \qquad (8) $$

where $\rho_{f0} = 1022$ kg m$^{-3}$ and $\mu_{f0} = 1.17 \times 10^{-3}$ Pa s are the density and viscosity values at the top of the model, $\beta = 5.9 \times 10^{-4}$ °C$^{-1}$ is the fluid thermal expansion, $T_t$ is the temperature at the top (20 °C), and $T_v = 62.1$ °C, is an approximate temperature to fit fluid viscosity in Eq. (8) (see Wang
<table><caption>**Table 3.** Boundary conditions used in Magri et al. (2017) for the 3D results.</caption><thead><tr><th>Variable</th><th>Top</th><th>Bottom</th><th>West</th><th>East</th><th>South</th><th>North</th></tr></thead><tbody><tr><td>Pore pressure</td><td>p<sub>f</sub> = 1 × 10<sup>5</sup> Pa</td><td>∇p<sub>f</sub> = 0</td><td>∇p<sub>f</sub> = 0</td><td>∇p<sub>f</sub> = 0</td><td>∇p<sub>f</sub> = 0</td><td>∇p<sub>f</sub> = 0</td></tr><tr><td>Temperature</td><td>T = 20 °C</td><td>T = 170 °C</td><td>∇T = 0</td><td>∇T = 0</td><td>∇T = 0</td><td>∇T = 0</td></tr></tbody></table>

<table><caption>**Table 4.** Values of the physical parameters used for the 3D benchmark according to Magri et al. (2017). Porosity values are on a scale from 0 to 1. Permeabilities are isotropic.</caption><thead><tr><th>Enclosing rocks</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Porosity <i>n</i></td><td>[1]</td><td>0.001</td></tr><tr><td>Density <i>ρ</i><sub>s</sub></td><td>kg m<sup>-3</sup></td><td>3125</td></tr><tr><td>Specific solid heat capacity <i>c</i><sub>s</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>800.0</td></tr><tr><td>Solid conductivity <i>λ</i><sub>s</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>3.0</td></tr><tr><td>Isotropic permeability <i>k</i></td><td>m<sup>2</sup></td><td>1 × 10<sup>-18</sup></td></tr><tr><td colspan="3"><strong>Fault</strong></td></tr><tr><td>Porosity <i>n</i></td><td>[1]</td><td>0.5</td></tr><tr><td>Solid density <i>ρ</i><sub>s</sub></td><td>kg m<sup>-3</sup></td><td>3125</td></tr><tr><td>Specific solid heat capacity <i>c</i><sub>s</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>800.0</td></tr><tr><td>Solid conductivity <i>λ</i><sub>s</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>3.0</td></tr><tr><td>Isotropic permeability <i>k</i></td><td>m<sup>2</sup></td><td>1 × 10<sup>-14</sup></td></tr><tr><td colspan="3"><strong>Fluid</strong></td></tr><tr><td>Top fluid density <i>ρ</i><sub>f0</sub></td><td>kg m<sup>-3</sup></td><td>1022</td></tr><tr><td>Top fluid viscosity <i>μ</i><sub>f0</sub></td><td>Pa s</td><td>0.0017</td></tr><tr><td>Fluid conductivity <i>λ</i><sub>f</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>0.65</td></tr><tr><td>Specific fluid heat capacity <i>c</i><sub>f</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>3925</td></tr></tbody></table>

et al., 1987, and Magri et al., 2017). The initial condition for temperature consists in a linear gradient from a temperature of 170 °C at the bottom to 20 °C at the top, therefore the gradient is approximately 27 °C km⁻¹. Additionally, a sinusoidal disturbance ε is added within the fault to trigger convection, and to have the same initial conditions proposed by Magri et al. (2017), as follows:

$$ \epsilon(x, y, z, t = 0) = \sin(\pi z) \cos(\pi y) \quad (9) $$

For the benchmark results, we follow Magri et al. (2017), but also Guillou-Frottier et al. (2020), where the 2017 3D benchmark was also reproduced. In the Fig. 4, we present a comparison between the thermal perturbation ($T - T(t = 0)$) within the fault at the time of 4027.14 kyr obtained in this work, and the results obtained in Magri et al. (2017). For this result we use an isotropic permeability of $1 \times 10^{-14}$ m² in the fault. This value differs from the permeability originally reported ($1.019 \times 10^{-13}$ m²) in Magri et al. (2017), but this value was not correct (Fabien Magri, personal communication, 2019). This permeability value ($1 \times 10^{-14}$ m²) was the best to reproduce the convective pattern and the thermal perturbation ($T - T(t = 0)$), which is in the range ($-21$ to

<table><caption>**Table 5.** Values taken from Magnenet et al. (2014) et Vallier et al. (2019) based on the Soultz-sous-Forêts reservoir. The granite has a higher permeability value than sediments, because it represents an altered granite. This value was obtained after an inversion developed at Vallier et al. (2019), permeability here is isotropic. Porosity values are on a scale from 0 to 1.</caption><tbody><tr><th colspan="3">Upper sediments</th></tr><tr><td>Porosity <i>n</i></td><td>[1]</td><td>0.19</td></tr><tr><td>Isotropic permeability <i>k</i></td><td>m<sup>2</sup></td><td>6 × 10<sup>-17</sup></td></tr><tr><td>Solid density <i>ρ</i><sub>s</sub></td><td>kg m<sup>-3</sup></td><td>2390.0</td></tr><tr><td>Specific solid heat capacity <i>c</i><sub>s</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>800.0</td></tr><tr><td>Solid thermal conductivity <i>λ</i><sub>s</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>0.6</td></tr><tr><th colspan="3">Lower sediments</th></tr><tr><td>Porosity <i>n</i></td><td>[1]</td><td>0.12</td></tr><tr><td>Isotropic permeability <i>k</i></td><td>m<sup>2</sup></td><td>6 × 10<sup>-17</sup></td></tr><tr><td>Solid density <i>ρ</i><sub>s</sub></td><td>kg m<sup>-3</sup></td><td>2390.0</td></tr><tr><td>Specific solid heat capacity <i>c</i><sub>s</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>800.0</td></tr><tr><td>Solid thermal conductivity <i>λ</i><sub>s</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>1.7</td></tr><tr><th colspan="3">Granite basement</th></tr><tr><td>Porosity <i>n</i></td><td>[1]</td><td>0.1</td></tr><tr><td>Isotropic permeability <i>k</i></td><td>m<sup>2</sup></td><td>4 × 10<sup>-14</sup></td></tr><tr><td>Solid density <i>ρ</i><sub>s</sub></td><td>kg m<sup>-3</sup></td><td>2630.0</td></tr><tr><td>Specific solid heat capacity <i>c</i><sub>s</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>800.0</td></tr><tr><td>Solid thermal conductivity <i>λ</i><sub>s</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>1.7</td></tr><tr><th colspan="3">Fluid</th></tr><tr><td>Top fluid density <i>ρ</i><sub>f0</sub></td><td>kg m<sup>-3</sup></td><td>1070</td></tr><tr><td>Specific fluid heat capacity <i>c</i><sub>f</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>4184</td></tr><tr><td>Fluid thermal conductivity <i>λ</i><sub>f</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>0.65</td></tr><tr><td>Top viscosity <i>μ</i><sub>f0</sub></td><td>Pa s</td><td>0.0017</td></tr></tbody></table>

30) °C in our work, and ($-20.62$ to $28.94$) °C in Magri et al. (2017).

## 4 Vendenheim site modeling

Figure 5 presents a simplified model of the Vendenheim site, which describes the main geological units in tabular form as a proxy for the sedimentary basin. These units include the granite basement, lower and upper sediments, which are located at depths deduced from the GEORG platform (Fig. 5). Each unit is homogenized and described by physical rock properties known at Soultz-sous-Forêts. Individual faults are
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison between the thermal perturbation in our benchmark result (a) and the results from Magri et al. (2017) (b), at a time of $t = 317.097$ kyr. Thermal perturbation is defined as: $T - T_0$, where $T_0$ is the temperature at time $t = 0$. The location of this result is in the fault (the red body defined in Fig. 3).</figcaption></figure>

not described as we search for large-scale hydro-thermal behavior. Previous studies (Vallier et al., 2019) have shown that regional faults have little impact on the local hydro-thermal behavior. The physical properties are taken from Vallier et al. (2019) (the name of the layers correspond to this article) and Magnenet et al. (2014), they are reported in Table 5. This properties are based on the units of the Soultz-sous-Forêts reservoir. The dimensions of the model are 6, 8, and 8 km in the x, y, and z directions respectively, with 100 m per side hexahedral elements. The upper boundary represents the Earth's surface.

The physical laws governing fluid density $ρ_f$ and viscosity $μ$ are identical to those presented in Magri et al. (2017) (Eq. 8). In this study we use the reported values of $ρ_{f0}$ and $μ_{f0}$, as presented in Table 5, to develop our analytical framework.

## Numerical solution

The boundary conditions exhibit no mass flow at all boundaries, except for the top boundary where the atmospheric pressure value is imposed. For temperature, heat flow is zero at all boundaries except for the top boundary where a value of 15 °C, that represents the average air temperature at the surface, is imposed. And at the bottom boundary, the system is heated with a value of 423 °C. This boundary conditions are presented in Table 6.

The distribution map of temperature and profiles using the properties of Table 5 are shown in Fig. 6 at the time of 4.77 kyr. We choose a time when the convective patterns have already formed. We observe a transient convective regime

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Simplified model based in the GEORG platform (https://maps.geopotenziale.eu, last access: 8 October 2024) over the Vendenheim site (North Strasbourg, France) and in the 2D model proposed in Vallier et al. (2019). The coordinates of the A, B, C, D, and central points, that refer to the location of the profiles of Fig. 6, are in kilometers respectively: (1.5,2), (4.5,2), (4.5,6), (1.5,6), (3,4).</figcaption></figure>

in the granite and a conductive regime in the sediments. The Rayleigh number in the granite, upper and lower sediments are 536 and 1.07, and 0.0289 respectively, the critical Rayleigh number being $4π^2 = 39.47$. Therefore, convection is expected in the granite but not in the sediments. In the vertical temperature profiles we added the measured tem-
<table><caption>**Table 6.** Boundary conditions used for the 3D results in the homogeneous model of Strasbourg.</caption><thead><tr><th>Variable</th><th>Top</th><th>Bottom</th><th>West</th><th>East</th><th>South</th><th>North</th></tr></thead><tbody><tr><td>Pore pressure</td><td>p<sub>f</sub> = 101 325 Pa</td><td>∇p<sub>f</sub> = 0</td><td>∇p<sub>f</sub> = 0</td><td>∇p<sub>f</sub> = 0</td><td>∇p<sub>f</sub> = 0</td><td>∇p<sub>f</sub> = 0</td></tr><tr><td>Temperature</td><td>T = 15 °C</td><td>T = 423 °C</td><td>∇T = 0</td><td>∇T = 0</td><td>∇T = 0</td><td>∇T = 0</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Temperature distribution in 4.77 kyr, with a convective regime in the granite and a conductive regime in the sediments. The horizontal width of the convection cells measure between 780 and 880 m. (b) Profiles showing the temperature along vertical lines at the different points A, B, C, D, Central, introduced in Fig. 5. We show the measured temperature profile at Cronembourg area. The Rayleigh number for the granite lower and upper sediments are 536, 0.0280, and 1.07 respectively.</figcaption></figure>

perature profile at Cronembourg site (at a distance of approximately 9 km from the Vendenheim site), measured by the “Bureau de Recherches Geologiques et Minieres” (BRGM) (personal communication, 1982), which is 5.10 °C per 100 m, showing that there exist a good fit between the calculated profiles and this gradient. The profiles show that the transition between the conductive and convective regimes occurs at a depth of about 5 km, according to Schmittbuhl et al. (2021) the earthquakes in the site are globally at a depth of about 5 km. This means that there may be a relationship between the transition region between the conductive and convective regimes and the region of the area where the induced earthquakes were generated.

## 5 Conclusions

We solve 2D and 3D benchmarks to test our methods based in TH coupled system, using the moose-based open-source application “Golem”. We show that for the 2D benchmark where the thickness is 1500 m, a transition from conductive to free convective solution exist for a permeability of the order of $1 \times 10^{-16} m^2$ with a Darcy velocity at a maximum of the order of several decimeters per year.

We developed a simplified model of the Vendenheim site, based on the regional geological model Georg, we apply a TH coupled system in this model, using the rock and fluid physical properties from previous studies at Soultz-sous-Forêts but adjusting their depths according to the features of the local geological model. We conclude that given the observed parameters at the Vendenheim site, a transient hydrothermal convective regime is observed at the granite, while a conductive behaviour is observed in the sediments. The approximate depth of induced seismic events (5 km) in the region coincides with the depth of the transition between the convective and conductive zones. However, as stated in the Introduction, exploring this possible relationship is beyond the scope of this article. Therefore, this work serves to establish the initial conditions for a future induced-seismicity model in the Vendenheim area.

*Code availability.* To calculate our results, we use the open source code Golem (https://git.gfz-potsdam.de, last access: 8 October 2024; https://doi.org/10.5281/zenodo.999401, Jacquey and Cacec, 2017).
*Data availability.* In this work, we validate our methods with two benchmarks, whose data (rock properties, boundary conditions and initial conditions) are reported in the references (Magri et al., 2017; Guillou-Frottier et al., 2020). We also studied the observed parameters of the lithologies of the Soults-sous-Forêts area, which is very close to the Vendenheim area, these data are reported in the references (Vallier et al., 2019).

*Author contributions.* JAT developed and performed the simulations, analysed the results and wrote the manuscript. GH, VM, MC and GB helped JAT in the simulation. MC, GB and JS contributed to the revision of the manuscript, giving their suggestions. JS contributed to analyze the results.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

*Special issue statement.* This article is part of the special issue "European Geosciences Union General Assembly 2024, EGU Division Energy, Resources & Environment (ERE)". It is a result of the EGU General Assembly 2024, Vienna, Austria, 14–19 April 2024.

*Acknowledgements.* This work of the Interdisciplinary Thematic Institut Geosicences for the energy system transition, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, was supported by Idex Unistra (ANR-10-IDEX-0002), and by the by the SFRI-STRAT'US project (ANR-20-SFRI-0012). We also thank the European Union's Horizon 2020 research and innovation program (DT-Geo, grant no. 101058129) for the funding of this work.

We are also grateful for the kind suggestions of the two anonymous reviewers who helped to improve this article.

*Financial support.* This work was funded by the ITI 2021-2028 programme of the University of Strasbourg, the CNRS and Inserm as a part of the Idex Unistra project (grant no. ANR-10-IDEX-0002) and by the by the SFRI-STRAT'US project (grant no. ANR-20-SFRI-0012). It has also been funded by the European Union's Horizon 2020 research and innovation programme (DT-Geo, grant no. 101058129).

*Review statement.* This paper was edited by Viktor J. Bruckman and reviewed by two anonymous referees.

## References

Angiboust, S., Wolf, S., Burov, E., Agard, P., and Yamato, P.: Effect of fluid circulation on subduction interface tectonic processes: Insights from thermo-mechanical numerical modelling, Earth Planet. Sc. Lett., 357, 238–248, 2012.

Bataillé, A., Genthon, P., Rabinowicz, M., and Fritz, B.: Modeling the coupling between free and forced convection in a vertical permeable slot: Implications for the heat production of an Enhanced Geothermal System, Geothermics, 35, 654–682, 2006.

Bayer, P., Attard, G., Blum, P., and Menberg, K.: The geothermal potential of cities, Renewable and Sustainable Energy Reviews, 106, 17–30, 2019.

Blöcher, G., Cacace, M., Jacquey, A. B., Zang, A., Heidbach, O., Hofmann, H., Kluge, C., and Zimmermann, G.: Evaluating micro-seismic events triggered by reservoir operations at the geothermal site of Groß Schönebeck (Germany), Rock Mech. Rock Eng., 51, 3265–3279, 2018.

Cacace, M. and Jacquey, A. B.: Flexible parallel implicit modelling of coupled thermal–hydraulic–mechanical processes in fractured rocks, Solid Earth, 8, 921–941, https://doi.org/10.5194/se-8-921-2017, 2017.

Davaille, A.: Two-layer thermal convection in miscible viscous fluids, J. Fluid Mech., 379, 223–253, 1999.

Desaive, T. C.: Thermoconvection dans les milieux poreux: stabilite et dynamique non-lineaire, Ph.D. thesis, Universite de Liege, Belgium, 2003.

Detournay, C., Damjanac, B., Torres, M., Cundall, P., Ligocki, L., and Gil, I.: Heat advection and forced convection in a lattice code–Implementation and geothermal applications, Rock Mechanics Bulletin, 1, 100004, https://doi.org/10.1016/j.rockmb.2022.100004, 2022.

Donaldson, I.: The simulation of geothermal systems with a simple convective model, Geothermics, 2, 649–654, 1970.

Duwiquet, H., Arbaret, L., Guillou-Frottier, L., Heap, M. J., and Bellanger, M.: On the geothermal potential of crustal fault zones: a case study from the Pontgibaud area (French Massif Central, France), Geothermal Energy, 7, 1–29, 2019.

Frey, M., Bär, K., Stober, I., Reinecker, J., van der Vaart, J., and Sass, I.: Assessment of deep geothermal research and development in the Upper Rhine Graben, Geothermal Energy, 10, 18, https://doi.org/10.1186/s40517-022-00226-2, 2022.

Gaucher, E., Schoenball, M., Heidbach, O., Zang, A., Fokker, P. A., van Wees, J.-D., and Kohl, T.: Induced seismicity in geothermal reservoirs: A review of forecasting approaches, Renewable and Sustainable Energy Reviews, 52, 1473–1490, 2015.

Giudicelli, G., Lindsay, A., Harbour, L., Icenhour, C., Li, M., Hansel, J. E., German, P., Behne, P., Marin, O., Stogner, R. H., Miller, J. M., Schwen, D., Wang, Y., Munday, L., Schunert, S., Spencer, B. W., Yushu, D., Recuero, A., Prince, Z. M., Nezdyur, M., Hu, T., Miao, Y., Jung, Y. S., Matthews, C., Novak, A., Langley, B., Truster, T., Nobre, N., Alger, B., Andrš, D., Kong, F., Carlsen, R., Slaughter, A. E., Peterson, J. W., Gaston, D., and Permann, C.: 3.0 – MOOSE: Enabling massively parallel multiphysics simulations, SoftwareX, 26, 101690, https://doi.org/10.1016/j.softx.2024.101690, 2024.

Guillou-Frottier, L., Duwiquet, H., Launay, G., Taillefer, A., Roche, V., and Link, G.: On the morphology and amplitude of 2D and 3D thermal anomalies induced by buoyancy-driven flow
within and around fault zones, Solid Earth, 11, 1571–1595,
https://doi.org/10.5194/se-11-1571-2020, 2020.

Jacquey, A. B. and Cacace, M.: GOLEM, a MOOSE-based application v1.0, Zenodo [code], https://doi.org/10.5281/zenodo.999401, 2017.

Kwiatek, G., Saarno, T., Ader, T., Bluemle, F., Bohnhoff, M., Chendorain, M., Dresen, G., Heikkinen, P., Kukkonen, I., Leary, P., Leonhardt, M., Malin, P., Martínez-Garzón, P., Passmore, K., Passmore, P., Valenzuela, S., Wollin, C.: Controlling fluid-induced seismicity during a 6.1-km-deep geothermal stimulation in Finland, Sci. Adv., 5, eaav7224, https://doi.org/10.1126/sciadv.aav7224, 2019.

Lapwood, E.: Convection of a fluid in a porous medium, in: Mathematical Proceedings of the Cambridge Philosophical Society, 44, 508–521, Cambridge University Press, 1948.

Liu, X., Wang, Y., Li, S., and Wang, J.: Convection heat transfer of supercritical CO₂ in a single fracture in enhanced geothermal systems, Int. Commun. Heat Mass, 123, 105170, https://doi.org/10.1016/j.icheatmasstransfer.2021.105170, 2021.

Lopez, T., Antoine, R., Kerr, Y., Darrozes, J., Rabinowicz, M., Ramillien, G., Cazenave, A., and Genthon, P.: Subsurface hydrology of the Lake Chad Basin from convection modelling and observations, Remote Sensing and Water Resources, 55, 281–312, 2016.

Magnenet, V., Fond, C., Genter, A., and Schmittbuhl, J.: Two-dimensional THM modelling of the large scale natural hydrothermal circulation at Soultz-sous-Forêts, Geothermal Energy, 2, 1–21, 2014.

Magri, F., Cacace, M., Fischer, T., Kolditz, O., Wang, W., and Watanabe, N.: Thermal convection of viscous fluids in a faulted system: 3D benchmark for numerical codes, Enrgy. Proced., 125, 310–317, 2017.

Majer, E. L., Baria, R., Stark, M., Oates, S., Bommer, J., Smith, B., and Asanuma, H.: Induced seismicity associated with enhanced geothermal systems, Geothermics, 36, 185–222, 2007.

Murphy, H. D.: Convective instabilities in vertical fractures and faults, J. Geophys. Res.-Sol. Ea., 84, 6121–6130, 1979.

Nield, D.: Estimation of an effective Rayleigh number for convection in a vertically inhomogeneous porous medium or clear fluid, Int. J. Heat Fluid Fl., 15, 337–340, 1994.

Patterson, J. W., Driesner, T., and Matthai, S. K.: Self-organizing fluid convection patterns in an en echelon fault array, Geophys. Res. Lett., 45, 4799–4808, 2018.

Pola, M., Cacace, M., Fabbri, P., Piccinini, L., Zampieri, D., and Torresan, F.: Fault control on a thermal anomaly: Conceptual and numerical modeling of a low-temperature geothermal system in the Southern Alps foreland basin (NE Italy), J. Geophys. Res.-Sol. Ea., 125, e2019JB017394, https://doi.org/10.1029/2019JB017394, 2020.

Ranalli, G. and Rybach, L.: Heat flow, heat transfer and lithosphere rheology in geothermal areas: features and examples, J. Volcanol. Geoth. Res., 148, 3–19, 2005.

Schmittbuhl, J., Lambotte, S., Lengliné, O., Grunberg, M., Jund, H., Vergne, J., Cornet, F., Doubre, C., and Masson, F.: Induced and triggered seismicity below the city of Strasbourg, France from November 2019 to January 2021, Comptes Rendus Géoscience, 353, 561–584, 2021.

Schubert, G. and Straus, J. M.: Three-dimensional and multicellular steady and unsteady convection in fluid-saturated porous media at high Rayleigh numbers, J. Fluid Mech., 94, 25–38, 1979.

Schubert, G., Turcotte, D. L., and Olson, P.: Mantle convection in the Earth and planets, Cambridge University Press, ISBN-10 052135367X, ISBN-13 978-0521353670, 2001.

Simms, M. and Garven, G.: Thermal convection in faulted extensional sedimentary basins: theoretical results from finite-element modeling, Geofluids, 4, 109–130, 2004.

Vallier, B., Magnenet, V., Schmittbuhl, J., and Fond, C.: Large scale hydro-thermal circulation in the deep geothermal reservoir of Soultz-sous-Forêts (France), Geothermics, 78, 154–169, 2019.

Wang, M., Kassoy, D., and Weidman, P.: Onset of convection in a vertical slab of saturated porous media between two impermeable conducting blocks, Int. J. Heat Mass Tran., 30, 1331–1341, 1987.

Zang, A., Oye, V., Jousset, P., Deichmann, N., Gritto, R., McGarr, A., Majer, E., and Bruhn, D.: Analysis of induced seismicity in geothermal reservoirs–An overview, Geothermics, 52, 6–21, 2014.