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

# Pressure propagation during hydraulic stimulation: case study of the 2000 stimulation at Soultz-sous-Forêts

Dariush Javani<fnref n="1" />* ID, Jean Schmittbuhl<fnref n="1" /> and François Cornet<fnref n="1" />

## Abstract

Hydraulic stimulation of pre-existing fractures and faults plays a significant role in improving hydraulic conductivity of the fracture network around injection and production wells in deep geothermal reservoirs. In present work, a three-dimensional distinct element method (3DEC, Itasca) is used to simulate the year 2000 hydraulic stimulation of GPK2 well of Soultz-sous-Forêts geothermal reservoir, where several major hydraulic stimulations have been performed and are well documented. The field scale numerical model of the reservoir (about 6000 × 4500 × 4500 m³) includes an explicit description of the main fault (FZ4770), was developed to constrain the large-scale hydromechanical properties of the fault, in particular, its behavior in terms of non-linear elastic response related to fault aperture changes. The first phase of the stimulation is modelled as a constant flow rate of 30 ls⁻<fnref n="1" /> of water injection into the center of a deformable fault at the depth of approximately 4.7 km. We observed that the fluid pressure front migration from the injection point along the fracture follows, under the in-situ stress condition and the moderate injection pressure, a pseudo-diffusion behavior as power-law function of time with a 0.5 exponent (vt) when the injection flow rate is constant. It is demonstrated that the dynamic evolution of aperture opening due to fluid injection into the fracture is responsible for the pressure propagation behavior, owing to a hydraulic aperture change rather than a fluid pressure diffusion process. This numerically observed propagation process is compatible with a high fault effective diffusivity of 13 m²/s as that observed in the field. In case of a linear increase of the injection flow rate, the pseudo-diffusion process disappears leading to a time dependent power-law migration of the pressure front with exponent of 0.75. The pressure propagation is shown to be strongly influenced by the injection scheme.

**Keywords:** Hydraulic stimulation, Deep fractured reservoirs, Soultz-sous-Forêts geothermal site, EGS

## Introduction

Hydraulic stimulation of pre-existing fractures and faults plays a significant role in improving hydraulic conductivity of the fracture network around injection and production wells in geothermal reservoirs. This operation is normally accompanied with

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

<footnote>*Correspondence:
javani.dariush@gmail.com</footnote>

<footnote n="1">¹ ITES/EOST, University
of Strasbourg/CNRS, Strasbourg,
France</footnote>
seismicity around the stimulation area which is classically supposed to be caused by fluid pressure diffusion through the fractured medium (William and Ellsworth 2013; Shapiro 2015; Cornet 2016; Orlecka-Sikora et al. 2020; Cauchie et al. 2020). Diffusion of fluid pressure, in consequence, may activate critically stressed fractures (or discontinuities close to slip) under in-situ stress state (Zoback 2007; Cornet et al. 2007). Shapiro et al. 2002, called this situation as near-failure equilibrium of rocks. In the other words, fluid pressure decreases fracture's inherent resistance against slip (Byerlee, 1978) by reducing effective normal stress acting on the fracture plane and causes fracture slip (Rice 1992; Shapiro et al. 1999; Parotidis et al. 2004; Barth et al. 2013; Cornet 2015; Blöcher et al. 2018); however, other mechanisms responsible for induced seismicity exist. Cappa et al. 2019, modelled a stepwise fluid injection in a fault and found that the increase in fluid pressure in the fault causes aseismic slip and fracture opening which accelerates by increase of the fluid pressure; indicating a stable fault creep during fluid injection. Dublanchet 2019, showed a two-phase response of the fault to a constant rate local fluid injection. He found that in the first step, the fluid pressure enhancement causes a reduction in effective normal stress which, in turn, triggers a slip solely in the pressurized zone and crack expands diffusively along strike. The shear crack in the second phase expands much faster than the fluid pressure diffusion; leading to, upon to the magnitude of initial shear stress along the fault compared to the residual shear strength left within the crack, an unstable (accelerating) crack regime or a stable constant rate crack expansion (which is the most probable case for tectonic faults).

Many authors (e.g. Ohtake 1974; Fletcher and Sykes 1977; Talwani and Acree 1985; Simpson et al. 1988) have estimated crust permeability or hydraulic diffusivity using seismicity-based methods. Shapiro et al. 1997 used fluid-injection-induced seismic emission at the KTB site, Germany, to estimate crust permeability. According to Shapiro et al. 1997, the physical background behind these estimations is the assumption that the stress state in the crust is close to a critical one such that the crust is in near-failure state. Hence, any small magnitude of stress perturbation, which can be generated by injected fluid, can disturb this equilibrium by decreasing effective normal stress. Thus, spatio-temporal onset of the seismicity relative to the start of the injection is interpreted as the time delay, necessary for fluid pressure to diffuse and generate a sufficiently large pressure perturbation in a distance to trigger seismicity. It is, therefore, hydraulic diffusivity that controls the spatio-temporal evolution of the microseismic events (Segal and Lu, 2015; Jin and Zoback 2017).

Linear diffusion of pore pressure in porous media has hydraulic diffusivity as a central parameter (Jaeger et al. 2007):

$$ \frac{\partial p}{\partial t} = D\nabla^2 p \qquad (1) $$

in which temporal evolution of the pore pressure ($\partial p/\partial t$) is related to the spatial evolution of the pore pressure ($\nabla^2 p$) by a proportionality factor of hydraulic diffusivity (D) (Jaeger et al. 2007; Wang 2000; Rozhko 2010). In fact, hydraulic diffusivity is a ratio of transport (permeability) to storage (specific storage capacity) properties, which in turn, depend on properties of rock (geometrical characteristics of flow conduits and deformation characteristics) and fluid (compressibility and viscosity) (Deng et al. 2021).
Shapiro et al. 1997, evaluated spatio-temporal distribution of approximately over 400 microseismic events induced during hydraulic fracturing experiment of KTB borehole in Germany to estimate permeability. At their study, fluid injection in KTB was approximated by a point source fluid injection into an infinite homogeneous, isotropic, poroelastic, and saturated medium for which the pore pressure diffusion can be considered as the mechanics of poroelastic media. Shapiro et al. 1997, solved pore pressure diffusion equation [Eq. (1)] and found that one solution of that can be as

$$r = \sqrt{4\pi Dt} \qquad (2)$$

where r is distance from injection point, D is hydraulic diffusivity of the medium, and t is time.

Murphy et al. 2004, showed that linear diffusion of fluid pressure is one possible scenario among three that may happen in a hydraulic stimulation only when injection rate and pressure are low such that joint aperture is unaffected. In the contrary of this condition, there is another scenario in which the joint aperture is evolving during injection and consequently to cause an effective hydraulic aperture change of the fault assuming a parallel plate model of the fault aperture but more realistically a change in hydraulic aperture due to the impact of the fluid injection on the sealing material along the fault in between the fault asperities that hold the normal load on the fault.

The aim of the present work is to numerically simulate two dimensional steady and laminar flow of a slightly compressible fluid in a fault modelled as two deformable parallel and smooth plates. The model is applied to a case study: the year 2000 hydraulic stimulation carried out on the well GPK2 of Soultz-sous-Forêts geothermal reservoir so as to characterize fluid pressure propagation through a large fault.

The Soultz-sous-Forêts geothermal reservoir is one of the first Enhanced Geothermal System (EGS) sites in the world (Genter et al., 2010). It hosts four operational wells, where GPK2, GPK3, and GPK4 being extended to a depth of 5000 m (Sausse et al. 2010). In the year 2000, the well GPK2 underwent a hydraulic stimulation to increase the hydraulic conductivity of the fracture system in the reservoir. It was a 6-day operation to inject 23,400 m3 water into open-hole section of well at depth of 4400 to 5050 m (Weidler et al. 2002). This injection procedure was accompanied by more than 24,000 microseismic events detected by surface and borehole receiver networks, of which 7200 were selected by Calo et al. (2011) for 4D P-wave tomography, where a variation in P-wave velocity was observed in the distance of over 1000 m in a period of about 3 h. Hydraulic data for this, as shown in Fig. 1, were obtained from the Data Center for Deep Geothermal Energy in Alsace (https://cdgp.u-strasbg.fr/ and https://tcs.ah-epos.eu/). Injection was started with 30 $ls^{-1}$ fluid injection rate for 24 h which is aimed to be simulated in the current study. This injection procedure is shown with blue curve in Fig. 1, where a fluid with constant flow rate of 30 $ls^{-1}$ is injected for 24 h; however, hydraulic stimulation was continued in two more steps with 40 and 50 $ls^{-1}$ for 27 and 90 h, respectively.

There is a well-known fracture at Soultz-sous-Forêts reservoir which intersects three wells constructing the EGS pilot plant (i.e., GPK2, GPK3, and GPK4). This fracture, shown in Fig. 3, is the largest fault in the reservoir and is designated GPK3-FZ4770 by Sausse et al. (2010). It has a dip direction and dip of 234° and 71° (64° in Dezayes et al. 2010) (Sausse et al. 2010), respectively. Here, we modeled this fault as a finite circular
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Stimulation procedure of well GPK2; red curve is field data (https://cdgp.u-strasbg.fr/), blue is numerically simulated one</figcaption></figure>

discontinuity with diameter of 3 km at the center of which the fluid is injected following
the procedure illustrated in Fig. 1.

**Numerical scheme**

**Numerical method**

The 3DEC code (Itasca CG, Inc.) uses the Discrete Element Method (DEM) introduced by Cundall (1988). This numerical technique is used to simulate in discontinuous media consisting of discontinuities and the intact rock blocks between them. Each block in 3DEC is divided into internal tetrahedral zones and nodes (Itasca 2016). The method uses an explicit time-marching procedure, where at each timestep, Newton's laws of motion are applied to all nodes to update the velocity, displacement and nodal force.

In 3DEC, discontinuities between the blocks are treated as boundary for which the forces are updated based on applied constitutive law. A contact detection algorithm is also deployed, and the integration of the law of motion is used to obtain new positions of blocks and contacts at the next timestep. 3DEC also has the ability to analyze fluid flow in fractures in a fully coupled hydromechanical analysis (Itasca 2016; Yin et al. 2020 and references therein), where the flow plane, or the contact between blocks, is a two-dimensional surface in which flow occurs through its thickness. The flow plane is assumed to be fully saturated and is discretized into triangular flow zones defined by three vertices or flow-knots (Itasca 2016). We assume in the present study that the permeability of the matrix is significantly smaller than that along the fracture and neglect possible leak-off effects.

**Coupled hydro-mechanical model**

In the present study, a linear relationship is assumed between aperture and normal effective stress, leading to a constant normal joint stiffness, which is the default behaviour in 3DEC,
as shown in Fig. 2. It is a simplification of the largely observed non-linear evolution of the closure with the normal stress (Bandis et al. 1983; Barton et al. 1985; Zangerl et al. 2008) for which a semi-logarithmic law is a good approximation for non-mated facing fracture surfaces. Such a semi-logarithmic law was shown to be related to a linear dependence of the stiffness with the normal stress or an exponential evolution of the stiffness with the fracture closure (Cook 1992). Recently, Deng et al. 2024 showed that such a behaviour was owing to the self-affine feature of fracture surfaces. For mated surfaces, a hyperbolic law is shown to better describe the closure behaviour (Bandis et al. 1983). The formulation for the latter case was tested in the PorePy framework (Ucar et al. 2017; Stefansson et al., 2024). Since a complex geometry of fracture surfaces is not implemented in the present model, we consider that a constant normal joint stiffness is a relevant approximation but further development will be addressed to account for such multi-scale geometry.

The constant normal stiffness of a fracture ($K_n$) means the existence of a linear relation between effective normal stress ($\sigma_{n'}$) and joint displacement ($d\delta$) (Hart 1993) as

$$ \sigma'_{n} = K_{n}d\delta \qquad (3) $$

The relation between hydraulic aperture and effective normal stress is a piecewise linear curve shown in Fig. 2. Hydraulic aperture ($h_H$) is obtained, in general, as

$$ h_H = h_{H0} + d\delta \qquad (4) $$

where $h_{H0}$ is the joint aperture at zero normal stress or initial hydraulic aperture.

The effect of changes in the fluid pressure of the fracture or in the applied stress yields using the concept of effective stress. According to Walsh (1981), the effect of pore pressure and applied normal stress on normal deformation can be expressed as

$$ \sigma'_{n} = \sigma_{n} - \beta p_{p} \qquad (5) $$

in which $p_p$ is pore pressure and $\beta$ is effective stress coefficient or Biot coefficient. Based on measurements of flow through simulated fractures carried out by Kranz et al. (1979),

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Modeled relationship between hydraulic aperture and effective normal stress for a rock joint (reproduced after Itasca 2016)</figcaption></figure>
$\beta$ can vary between 0.5 and 1.0 for tensile fractures with rough surfaces and fractures with smooth surfaces, respectively. For smooth fracture surfaces of current study, $\beta$ was considered to be unity, so that, Biot's effective stress becomes equivalent to Terzaghi's effective stress (Cornet 2015).

Updated hydraulic aperture values by Eqs. (3) and (4) are then substituted in cubic law (Witherspoon et al. 1980):

$$q = -\frac{h_H^3 \cdot w}{12 \mu_d} p_i \qquad (6)$$

To obtain the fluid flow rate through the contact or flow plane (q). In Eq. (6), w is flow plane width, $\mu_d$ is dynamic viscosity of fluid and $p_i$ is the pressure gradient. As mentioned earlier, the flow plane consists of triangular flow zones and three vertices or flow-knots, where the fluid pressure is stored and can communicate with neighboring flow-knots. The fluid pressure difference in the adjacent flow-knots is the driving mechanism of fluid flow along the flow zones. At each timestep of ($\Delta t$), the sum of fluid flow entering a flow-knot (Q) from all surrounding contacts causes a volume change of ($\Delta V$) and a new flow-knot pressure of p as below

$$p = p_0 + \frac{Q \Delta t}{c_f V} - \frac{1}{c_f} \frac{\Delta V}{V_m} \qquad (7)$$

where $p_0$ is the pressure value in the preceding time-step, $c_f$ is the fluid compressibility, and

$$\Delta V = V - V_0 \qquad (8)$$

$$V_m = (V + V_0)/2 \qquad (9)$$

where V and $V_0$ are flow-knot volume in two successive timesteps.

## Large-scale hydro-mechanical model

This section aims to describe the simulation of the first stage of hydraulic stimulation in Fig. 1 (blue line), where brine (fluid properties are summarized in Table 2) was injected with constant flow rate of 30 $\text{ls}^{-1}$ for 24 h.

### Model setup

Figure 3 illustrates the numerical model including a block of granite rock matrix containing a 3000 m diameter circular fracture called GPK3-FZ4770 with dip 71°, dip direction of 234° (Sausse et al. 2010) as well as the injection point located approximately at the center of the fracture at depth of 4.7 km. The fracture cuts the block in two parts as shown by light blue and light green in Fig. 3. These two parts have identical elastic, isotropic, and homogeneous properties.

In-situ stress state is determined by Valley and Evans, 2007, and Cornet et al. 2007, based on the data acquired from breakouts and drilling induced tensile fractures (DITFs) in the former and hydraulic stimulation of the year 1993 of GPK1 well. Pore pressure in the reservoir has also been considered to be hydrostatic in both references.
<figure><img src="image_5.png" /><figcaption>**Fig. 3** a Sketch of the 3DEC model around the GPK2 well with the embedded FZ4770 fracture (red circle with diameter of 3 km) and the injection point at 4700 m depth; b orientation of the in-situ stresses and fault</figcaption></figure>

<caption>**Table 1** In-situ stress state at Soultz-sous-Forêts; Z is depth in "m", stresses are in MPa (Cornet et al. 2007)</caption>

**In-situ stress state**

$$S_v = 33.8 + 0.0255(Z - 1377)$$

$$S_{hmin}/S_v = 0.54$$

$$0.95 < S_{hmax}S_v < 1.1$$

$$P_p = 0.9 + 0.0098Z$$

In present work, in-situ stress state introduced by Cornet et al. 2007, was taken into account as shown in Table 1, where $S_v$, $S_{hmin}$, and $S_{hmax}$ are, respectively, vertical, minimum horizontal, and maximum horizontal in-situ stresses in MPa, $P_p$ is pore pressure (MPa) and Z is depth in meters. According to Cornet et al. 2007, the maximum horizontal stress has a direction of N170°E ± 10° (169° ± 14° based on Klee and Rummel, 1993, and Valley and Evans, 2007). The relation between the vertical stress and the maximum horizontal stress in Table 1 indicates that the stress regime in the Soultz-sous-Forêts is normal to strike-slip according to the faulting regime classification of Anderson (1905). Here, it is assumed that the faulting regime is normal with a ratio of maximum horizontal stress to vertical stress of 0.95. Physical and mechanical properties of the rock matrix along the open-hole section of the well and the fracture properties were derived from existing literature and summarized in Table 2.

For elastic moduli of the rock matrix at Soultz-sous-Forêts reservoir, Vallier et al. (2019) did an extensive inversion of observed temperature and stress-depth profiles using a thermo-hydro-mechanical (THM) model to obtain the possible magnitudes of Poisson's ratio, and Young's modulus. They concluded that 25 GPa is the more realistic value for Young's modulus at large scale (100 m) of the lower granite, where the open-hole section of GPK2 is located. The density of the lower granite was also chosen to be 2680 kg/m³, according to Gentier et al. (2005). Gentier et al. (2005) carried out a numerical study using 3DEC on hydraulic stimulation of wells at Soultz-sous-Forêts to reproduce observed behavior in terms of pressure-flow curves and location of the stimulated fractures and found a good agreement between obtained numerical fluid flow results
<table><caption>**Table 2** Properties of rock, fracture, and fluid</caption><thead><tr><th>Element</th><th>Property</th><th>Symbol</th><th>Value</th><th>Unit</th><th>Remark</th></tr></thead><tbody><tr><td rowspan="3">Rock Matrix</td><td>Young's modulus</td><td>E</td><td>25*–47<br>52<br>25*</td><td>GPa</td><td>Rummel (1991)<br>Gentier et al. (2005)<br>Sausse (2002), Vallier et al. (2019)</td></tr><tr><td>Poisson's ratio</td><td>u</td><td>0.29*<br>0.1<br>0.38<br>0.20</td><td>–</td><td>Gentier et al. (2003) &amp; (2005)<br>Sausse (2002)<br>Rummel (1991)<br>Vallier et al. (2019)</td></tr><tr><td>Bulk density</td><td>ρ</td><td>2680*<br>2650–2800</td><td>kg/m³</td><td>Gentier et al. (2005)<br>Haenel (1983)</td></tr><tr><td rowspan="5">Fracture (inside area)</td><td>Shear stiffness</td><td>K<sub>s</sub></td><td>20*<br>80–180</td><td>GPa/m</td><td>Gentier et al. (2005)<br>Gentier et al. (2003)</td></tr><tr><td>Normal stiffness</td><td>K<sub>n</sub></td><td>80*<br>80–180</td><td>GPa/m</td><td>Gentier et al. (2005)<br>Gentier et al. (2003)</td></tr><tr><td>Friction angle</td><td>φ</td><td>39*<br>45<br>48</td><td>°</td><td>Cornet et al. (2007)<br>Evans (2005)<br>Gentier et al. (2005)</td></tr><tr><td>Cohesion</td><td>S<sub>0</sub></td><td>4–8<br>6.5*<br>0</td><td>MPa</td><td>Cornet et al. (2007)<br>Evans (2005)<br>Gentier et al. (2005)</td></tr><tr><td>Initial aperture</td><td>h<sub>0</sub></td><td>1.4e–04*<br>6e–07</td><td>m</td><td>Gentier et al. (2005)<br>Gentier et al. (2003)</td></tr><tr><td rowspan="3">Fracture (outside area)</td><td>Residual aperture</td><td>h<sub>res</sub></td><td>0.7e–04*<br>1e–08</td><td>m</td><td>Gentier et al. (2005)<br>Gentier et al. (2003)</td></tr><tr><td>Maximum aperture</td><td>h<sub>max</sub></td><td>120e–04*<br>1e–03</td><td>m</td><td>Gentier et al. (2005)<br>Gentier et al. (2003)</td></tr><tr><td>Shear stiffness</td><td>K<sub>s</sub></td><td>20*</td><td>GPa/m</td><td>Gentier et al. (2005)</td></tr><tr><td rowspan="5">Fluid</td><td>Normal stiffness</td><td>K<sub>n</sub></td><td>80*</td><td>GPa/m</td><td>Gentier et al. (2005)</td></tr><tr><td>Friction angle</td><td>φ</td><td>39*</td><td>°</td><td>Cornet et al. (2007)</td></tr><tr><td>Cohesion</td><td>S<sub>0</sub></td><td>6.5*</td><td>MPa</td><td>Cornet et al. (2007)</td></tr><tr><td>Initial aperture</td><td>h<sub>0</sub></td><td>1.4e–04*<br>1.4e–04*</td><td>m</td><td>Gentier et al. (2005)<br>Gentier et al. (2005)</td></tr><tr><td>Residual aperture</td><td>h<sub>res</sub></td><td>1.4e–04*</td><td>m</td><td>Gentier et al. (2005)</td></tr><tr><td></td><td>Maximum aperture</td><td>h<sub>max</sub></td><td>1.4e–04*</td><td>m</td><td>Gentier et al. (2005)</td></tr><tr><td></td><td>Density</td><td>ρ<sub>f</sub></td><td>1000*</td><td>kg/m³</td><td>Ball (2008)</td></tr><tr><td></td><td>Viscosity</td><td>μ<sub>f</sub></td><td>1*</td><td>Centpoise</td><td>Kestin et al. (1978)</td></tr><tr><td></td><td>Compressibility</td><td>c<sub>f</sub></td><td>5e–10*</td><td>1/Pa</td><td>Fine and Millero (1973)</td></tr></tbody></table>

The * indicates the values used in the current study

and the data achieved from in-situ hydraulic tests. The magnitudes of Poisson's ratio, fracture's normal and shear stiffnesses, initial, residual, and maximum apertures of the fracture which are those used by Gentier et al. (2005) in their models (Table 2). To have an appropriate value of fracture's cohesion and friction angle among those reported in the literature, we used the data provided by Cornet et al. (2007) in which the authors extensively discussed two hypotheses about shear failure of the pre-existing fractures and initiation of fresh fracture surfaces at Soultz-sous-Forêts reservoir. They stated that for pre-existing fractures to be just in equilibrium at the present stress state, or in other words, to be stimulated with a slight increase in in-situ pore pressure, the cohesion and friction angle of fractures should be, respectively, equal to 6.5 MPa and 39° (i.e., a friction coefficient of 0.81).

The model is discretized with tetrahedrons of side-size 200 m (mesh size of 200 m). In terms of boundary conditions, the bottom of the model is constrained by prescribing a
zero velocity boundary condition. For the sides of the model, a boundary stress condi-
tion to reproduce the in-situ stress state (Table 1) is applied. The model is first allowed
to reach mechanical equilibrium, then the injection procedure is applied to the center
of the fracture (assuming the fracture intersection with the wellbore is the center of the
fracture) and the hydromechanical calculation is performed.

Results

Fluid pressure evolution

By simulating the injection procedure illustrated by the blue line in Fig. 1, we observed, as shown in Fig. 4, that the fluid pressure at injection point raised sharply to a fluid pressure of approximately 13 MPa and after an approximately plateau phase steadily increased to a maximum value of over 20 MPa at the end of injection; which is somehow larger than observed field pressure during the first step of GPK2 stimulation at year 2000 (≈12 MPa), even though the sharp rise of pressure and plateau phase are similar to that of field observations. Numerical fluid pressure shown on y-axis of Fig. 4 is a proxy of well-head pressure (WHP) by considering the fact that in present effects on fluid pressure such as viscous drag, turbulent effect, etc. are negligible. Therefore, a proxy of the WHP was obtained by subtracting the hydrostatic pressure at injection point from downhole pressure (DHP). This might be one reason behind the observed difference between the numerical fluid pressure and its field counterpart among others, such as initial hydraulic aperture, permeability of surrounding rock, etc. In addition, it has to be noted that, in the current study, there is no leak-off in the model (i.e., undrained boundary condition) which might contribute to a relative over-pressure compared to the observed field pressure.

To investigate the fluid pressure propagation through the fracture, a set of measur-
ing points along a vertical cross section of the fault from the top of the fault (point

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Numerical WHP evolution in time at injection point in comparison with that observed in the field</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Locations of the measuring points (1:15) on the fracture, where the sphere number 8 is injection point; the triangles are flow zone and corresponding flow-knots</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Illustration of pressure build-up onset time determination at measuring point 4</figcaption></figure>

1) to the bottom of the fault (point 15) passing through the injection point (point 8), as shown in Fig. 5. These measuring points were used to pick a time (called pressure build-up onset time (PBOT), afterwards) at which the fluid pressure reached a pre-defined value which we set it as in-situ fluid pressure plus 1 MPa. Figure 6 illustrates fluid pressure evolution at a measuring point 4 located 886 m away from the injection point at which overpressure of 1 MPa on the y-axis is opted to find its corresponding PBOT. These PBOTs at various measuring points were then plotted against depths of the corresponding measuring points, as shown in Fig. 7, which is a plot of depth versus PBOT on the left-hand side (LHS) and distance from injection point versus PBOT on the right-hand side (RHS). It can be seen on the RHS graphs of Fig. 8 that the upward and downward velocities of the fluid pressure propagation are approximately equal and obeying a power-law function of time with an exponent of 0.5 as Eqs. (10) and (11) for, respectively, upward and downward propagation of fluid pressure:
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Left: pressure propagation from injection point to the top and bottom of the fracture on which the red curves are pressure evolution at each measuring point; Right: PBOT with respect to the distance of measuring points from the injection point to the shallower depth (top graph) and to the dipper depth (bottom graph) relative to the injection point</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8** Superimposition of WHP curves vs time at measuring points starting from overpressure of 1 MPa</figcaption></figure>
$$d_{p(t)} = 13.16t^{0.50} \qquad (10)$$

$$d_{p(t)} = 12.18t^{0.50} \qquad (11)$$

where $d_{p(t)}$ is pressure build-up distance in m and $t$ is time in second (these units are effective for Eqs. 12–17 as well). The red curves in the LHS graph in Fig. 7 represent the fluid pressure evolution at each measuring point starting from PBOT. Superimposition of these red curves of Fig. 7 are shown in Fig. 8 in which the time period starts from the PBOT at measuring points as determined using the process described in Fig. 6; so that, all the pressure curves start from 1 MPa.

Comparison of Eqs. (10) and (11) with Eq. (2) reveals that the coefficients of 13.16 or 12.18 in Eqs. (10) and (11) are a proxy of $\sqrt{4\pi D}$ of Eq. (2) which leads to an average effective hydraulic diffusivity of $D \sim 13 \text{ m}^2/\text{s}$. A closer look at the right sub-figures of Fig. 7 also shows a fluid pressure propagation of over 1000 m in 3 h (or 10,500 s), consistent with the observation of Calo et al. (2011), where a rapid migration of over 1000 m in 3 h was observed in the P-wave velocity variation during injection (Fig. 9).

<figure><img src="image_11.png" /><figcaption>**Fig. 9** Left: propagation of hydraulic aperture increase from injection point to the top and bottom of the fracture on which the red curves are hydraulic aperture evolution at each measuring point; Right: hydraulic aperture increase with respect to the distance of measuring points from the injection point to the shallower depth (top graph) and to the dipper depth (bottom graph) relative to the injection point</figcaption></figure>
## Hydraulic aperture evolution

Figure 9 illustrates velocity of fracture opening migration along the fracture towards upper and lower edges along its diameter. The onset of aperture increase was considered as the time at which the amount of the increase at hydraulic aperture exceeded the 5% of the initial value of that. Here, knowing the initial hydraulic aperture of 1.4e-04 m, the onset time was the moment at which the hydraulic aperture exceeded 1.47e-04 m. This process is demonstrated in Fig. 10. To avoid any boundary condition effect on hydraulic aperture variation (such as accumulation of the injected fluid at the edge of the fracture which causes early opening increase at the measuring points close to the fracture edge compared to those far from the edge), the data of the two last measuring points at each side of the fracture diameter were not considered for evaluations of hydraulic aperture variation in time. As it can be seen in Fig. 10 (sub-figures at the RHS), the propagation velocity of increase in the hydraulic aperture obeys a power-law function with an exponent of 0.5, similar to that of the fluid pressure propagation as

$$d_{a(t)} = 15.76t^{0.49} \qquad (12)$$

$$d_{a(t)} = 15.30t^{0.50} \qquad (13)$$

where $d_{a(t)}$ is the fracture opening distance relative to injection point.

To better characterize the response of fracture to fluid injection, a time-dependent flow rate of injection as $Q(t) = 1 \times 10^{-5} \times t$, in which $Q(t)$ is linearly increasing flow rate with time ($t$), was applied to the injection point so as to assess the effect of dynamic injection process on the propagation of fluid pressure as well as the fracture opening along the fracture. In this complementary test, the flow rate increased by $1 \times 10^{-5} \text{ m}^3$ at each second compared to the fixed flow rate of 0.03 $\text{m}^3/\text{s}$ in the reference case. A fluid volume of 2592 $\text{m}^3$, as equal as that of the constant injection rate case, was injected into the fracture for 22,768 s. Fluid pressure propagation profile and corresponding propagation velocity profiles were then plotted as shown in Fig. 11 in which a power-law relation between pressure propagation distance, m, and time, s, are obvious as

<figure><img src="image_12.png" /><figcaption>Fig. 10 Illustration of opening onset time selection at measuring point 4 shown in Fig. 6</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 11** Left: pressure propagation from injection point to the top and bottom of the fracture (time-dependent injection flow rate case); Right: PBOT with respect to the distance of measuring points from the injection point to the shallower depth (top graph) and to the dipper depth (bottom graph) relative to the injection point</figcaption></figure>

$$d_{p(t)} = 1.41t^{0.75} \qquad (14)$$

and

$$d_{p(t)} = 1.38t^{0.75} \qquad (15)$$

for, respectively, upward and downward profiles. An approximately similar trend is also exhibited by hydraulic aperture opening propagation along the fault plane as illustrated in Fig. 12 as

$$d_{a(t)} = 2.21t^{0.71} \qquad (16)$$

and

$$d_{a(t)} = 1.99t^{0.72} \qquad (17)$$

upward and downward, respectively.

According to these results, there exists a pressure front as well as a hydraulic aperture front in the constant flow rate simulation which both propagated in a square root diffusion-like trend. Existence of observed fronts in the fluid pressure and hydraulic aperture were
<figure><img src="image_14.png" /><figcaption>**Fig. 12** Left: propagation of hydraulic aperture increase from injection point to the top and bottom of the fracture (time-dependent injection flow rate case); Right: hydraulic aperture increase with respect to the distance of measuring points from the injection point to the shallower depth (top graph) and to the dipper depth (bottom graph) relative to the injection point</figcaption></figure>

verified through a complementary simulation, where the flow rate was linearly increased in time and, in consequence, propagation of fluid pressure and fracture opening along the fault turned into power-law function of time with exponent of 0.75.

## Discussion

### Comparison with a simple on-off opening model

Lengliné et al., 2017, proposed a non-diffusive model for the migration of the seismic activity occurring during the hydraulic stimulation at Ritthershoffen (close to the Soultz-sous-Forêts reservoir). This model assumes a piecewise Heaviside function for the fracture opening at the fracture front: the aperture is assumed to jump from zero (or $h_{H0}$) to $h_H$ with a very low stiffness effect once the pressure/overpressure exists in the fracture. The front position is supposed to move as, $r(t)$, separating a domain of an open fault with constant hydraulic aperture, $h_H$, and a closed domain with a zero aperture. This on-off opening model was expressed as

$$r(t) = \sqrt{\frac{\int_{0}^{t} Q(t')dt'}{\pi h_H}} \qquad (18)$$
where

$$V(t) = \int_{0}^{t} Q(t')dt' \qquad (19)$$

is the injected volume at evolving flow rates of $Q(t')$. Supposing that the injected fluid migrates over a disk of radius $r$, then injected volume can be calculated as

$$V(t) = \pi h_H r(t)^2 \qquad (20)$$

Substituting (20) in left hand side of (19) and reordering for $r(t)$ will lead to (18). Equation (18) predicts the radius of pressure propagation by assuming a constant aperture throughout the fracture in which the fluid flows. Application of this model to the observed seismic data of the Rittershoffen field resulted in an estimate of the seismic migration front that was at least as good a fit to the filed data as the classical diffusion interpretation of Shapiro et al. 1997.

For our case of study, it was attempted to test the on-off concept introduced by Lengliné et al., 2017 [Eq. (20)], using the results obtained from 3DEC simulations to have a prediction of pressure propagation for problems in which the fracture aperture is variable with pressure following stress equilibrium. As it can be seen in Eq. (20), the on-off model is a volumetric proxy for pressure propagation based on the injected fluid migration along the fracture, whereas 3DEC is a mechanically balanced solution based on hydromechanical coupling and effective stress equilibrium. Here, the point is to test if a volumetric approximation is a good proxy or not. For that, the average value of hydraulic aperture profile throughout the fracture at each time-step of numerical modeling, $h_{Ha}(t)$, was used, as a crude approximation of aperture in position and time, to determine the radius of pressurized zone.

Figure 13 illustrates hydraulic aperture profile along the fracture length as well as the average hydraulic aperture over the whole fracture plane for constant flow rate cases at 14,284 s after start of injection (it is the time when the pressure build-up reaches to the last measuring points). The average hydraulic aperture at each time step of the constant and time-varying injection flow rate cases is computed by a FISH code embedded in the developed numerical model. This average hydraulic aperture, calculated over the flow-knots located in the fracture surface, is used as a quasi-constant proxy for hydraulic aperture to be substituted into Eq. (18) for estimating the pressure migration distance versus time. The resulting estimates are shown in Figs. 14 and 15. Plots are made for the period when the overpressure reaches the furthest measurement point on the fracture plane.

As it can be seen, in both cases shown in Figs. 14 and 15, the radius predicted by the on-off opening model (Lengliné et al. 2017) is in good agreement with the

<figure><img src="image_15.png" /><figcaption>Fig. 13 Hydraulic aperture profile along the fracture at 14,284 s after start of injection (Solid line) and average hydraulic aperture (dashed line) over the fracture surface</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 14** Pressure front propagation vs time for constant flow rate case based on the "on-off opening model" (solid lines) and 3DEC model (dots)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15** Pressure front propagation vs time for variable flow rate case based on the "on-off opening model" (solid lines) and 3DEC model (dots)</figcaption></figure>

overpressure radius obtained with the present 3DEC model. This supports the hypothesis that in the studied hydraulic stimulation there is a fluid pressure front controlled by the hydraulic opening front, instead of being a linear pressure diffusion process. In addition, and as discussed earlier, the on-off opening model is a volumetric proxy for pressure migration along the fracture, whereas the 3DEC is a mechanically balanced solution. Looking at Figs. 14 and 15, where a good agreement between the 3DEC results and the lift-off curves is shown, it becomes clear that the volumetric approximation is a good proxy for characterizing the fluid pressure propagation through the fracture.
Simulation of the whole stimulation

As mentioned above, Calo et al 2011 reported a P-wave velocity evolution over a dis-
tance of more than 1000 m in a period of about 3 h during the hydraulic stimulation
of the GPK well in 2000. This is also shown by our numerical model results in the right
subfigures of Figs. 7 and 9 in which both the propagation of fluid pressure and fracture
opening fronts reached to a distance of over 1000 m at time of about 11,000 s (i.e., 3
h). These results were obtained by simulating the first phase of fluid injection in the
undrained boundary condition, where there was no leakage from the fracture wall and
the edge into the surrounding rock matrix. A series of complementary tests were per-
formed to simulate the full fluid injection process (shown in Fig. 16) to assess the per-
formance of the numerical models to cover the observed wellhead pressure evolution
during the injection phases and shut-in phase. This was done by first changing the und-
rained boundary condition to the drained condition, where the fluid was allowed to flow
through the outer plane up to the boundaries of the block, through which a discharge
equal to the injected flow rate was set to complete the drained condition in the model.
The parameters controlling the response of the fracture to the injected fluid were then
adjusted to mimic the observed pressure evolution in the field. Figure 17 shows the well-
head pressure evolution obtained from the numerical simulation compared to the field
observation. The numerical curve is obtained from the model, where the normal stiff-
ness of the fracture is set to 40 and 80 GPa/m inside and outside, respectively. In addi-
tion, the hydraulic aperture values in the outside of the fracture were as equal as those in
the inside of the fracture.

As can be seen from Fig. 17, the numerical wellhead pressure generally follows the
field observations, with an approximate 1 MPa sharp increase when the injection rate
is increased by 10 ls⁻¹, a continuous and gradual increase in the third stimulation stage,
and finally a 4 MPa sudden pressure drop in the shut-in phase and then a gradual decline
to a pressure of 5–6 MPa. However, the slope of the pressure drop in the shut-in phase is

<figure><img src="image_18.png" /><figcaption>**Fig. 16** Injection rate and pressure measured at GPK2 wellhead</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 17** Numerical WHP evolution over time at the injection point compared to that measured at the wellhead of the GPK2 well</figcaption></figure>

smaller in the field observation than in the numerical simulation, suggesting a long-term pressure decline in the field observation, which is addressed by Weidler et al. (2002) as a result of the lower rock mass permeability of the deeper reservoir compared to the upper reservoir of Soultz-sous-Forêts. In addition, and among many others, the more complex geometry of the fracture network compared to a single fracture may be another reason for the observed difference between the shut-in phase pressure drop in the field observation and the numerical simulation. An extended modeling of the post-injection phase including a single fault or a complex network of faults was tested by Ucar et al., (2017) who showed the importance of the non-constant normal stiffness during fracture closure. Furthermore, during the first and second stimulation stages of the field observation, a well-head pressure of more than 12 MPa was observed which was dropped for about 1 MPa (as can be seen in Figs. 16 and 17), but was not apparent in the numerical results, although the pressure curve rose sharply. This could be due to the fact that the fault behavior in the numerical simulation was elastic, while the field observation, as could be inferred from the recorded induced seismic events, occurred in the elastic-plastic regime. In fact, as explained by Weidler et al. (2002) about the field observation, in the first and second stages, after reaching the peak break-down pressure of about 12 MPa, the fluid flows along a more permeable medium and thus the pressure decreases continuously, while in the third stage the fluid pressure front has reached an impermeable outer boundary through which the fluid pressure rises from 12 to 14 MPa.

As mentioned above, the 2000 stimulation operation of GPK2 was accompanied by more than 24,000 microseismic events detected by surface and borehole receiver networks, of which 7200 were selected by Calo et al. (2011) for 4D P-wave tomography. The entire injection operation was divided into 14 sets based on the changes in flow rate and fluid pressure that occurred (Fig. 18), and 4D P-wave tomography was performed on the
<figure><img src="image_20.png" /><figcaption>**Fig. 18** Numerical WHP evolution over time at the injection point compared to that measured at the wellhead of the GPK2 well [reproduced after Calo et al. (2011)]</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 19** A Evolution of the P-wave seismic velocity at 4.6 km depth during the 2000 stimulation test (after Calo et al. 2011) on which the $V_p$ measuring point is shown by a yellow circle; B approximate P-wave velocity changes at the yellow point shown in A versus set and flow rate</figcaption></figure>

events that occurred in each set. Figure 19 shows 6 sets out of these 14, where a variation in P-wave velocity can be observed at a distance of over 1000 m over a period of approximately 3 h. A measurement point close to the injection source has been selected on the subplots of Fig. 19A, where an approximation of the P-wave velocity evolution has been plotted against the change in set and flow rate, as shown in Fig. 19B. It can be seen that an increase in P-wave velocity occurs as the flow rate is increased from 30 to 40 and from 40 to $50 \text{ ls}^{-1}$ when moving from set 2 to 3 and set 5 to 6, respectively. This behavior is inverse to the fact that any increase in pore/fluid pressure will decrease the P-wave velocity (Mavko and Vanorio 2010; Darot and Reuchlé, 2000; Todd and Simmons 1972). Considering the low permeability medium addressed by Weidler et al. (2002), in
addition to the existence of a fluid pressure/fracture opening front that propagates along the fracture faster than fluid diffusion in the fractured surrounding medium which was demonstrated in this study, one can state that the overpressure generated by the injected fluid in the fracture acts as an additive component to the confining pressure of the rock mass, instead of playing the role of pore pressure enhancement. In this situation, and according to the literature already mentioned, especially Todd and Simmons (1972), any increase in confining pressure can cause an increase in P-wave velocity. In our case, an overpressure of 1 MPa acting normal to the fracture plane with a dip and azimuth of 71° and 144°, respectively, can be translated into a confining pressure component of about 0.6 MPa on average in all three Cartesian directions, which, based on the results of Todd and Simmons, can produce a 2% increase in P-wave velocity, assuming an average confining pressure of 90 MPa and a pore pressure of about 50 MPa. This increase in P-wave velocity is comparable to the 10% and 2% increases seen in the field observations after the first and second injection rate increases, respectively, as can be roughly calculated from Fig. 19B.

## Conclusion

Numerical simulation of the large-scale hydraulic stimulation of the year 2000 performed at GPK2 well of the Soultz-sous-Forêts site is carried out to better characterize fluid pressure propagation. The 3DEC distinct element method of Itasca was used to simulate the behavior of the main fracture (FZ4770 fracture) of the reservoir (Sausse et al 2010). This main fracture of the reservoir is considered to be the dominant one and modeled as embedded in an elastic, homogeneous, isotropic, and impermeable block of granite during the injection of water for 24 h at a constant flow rate (30 $ls^{-1}$). The results indicated that a fluid pressure front migrates as a function of $\sqrt{t}$ when a constant injection rate is imposed with an average effective hydraulic diffusivity of ~13 $m^2/s$. This fluid pressure migrates a distance of over 1000 m in about 10,500 s consistent with the observations of Calo et al. 2011, where a rapid change in seismic velocity of over 1000 m in 3 h was observed; however, this value of hydraulic diffusivity is much larger than that estimated by Shapiro et al. 2002, based on the classical pore pressure diffusion process introduced by Shapiro et al. 1997 to describe seismicity migration. Interestingly, the same $\sqrt{t}$ behavior is observed in the fracture's hydraulic aperture. Application of a linearly increasing injection rate resulted in different behavior in both pressure and fracture aperture migration ($t^{0.75}$), confirming that the pressure front migration is not driven by a pressure diffusion process.

The findings of both aforementioned injection scenarios were compared to the model introduced by Lengliné et al., 2017. This model is a simple on-off opening model assuming a constant aperture along the open part of the crack and a sharp front between open and close part of the crack with no diffusion of the pressure. It expresses that it is the aperture evolution which is the driving mechanism of seismicity, and not the pressure diffusion. We generalized the model to our numerical results using the fracture aperture values computed from the 3DEC simulations during injection period. The outcome showed a good agreement between the two models indicating the fact that there exists a front in fluid pressure related to the hydraulic aperture front that propagates in a diffusion-like trend along the fracture plane. This rapidly propagating pressure front
along the fault, embedded in a low permeability medium, where fluid diffusion is slow compared to the propagation of the pressure front, can generate an additional confining pressure component that exerts more stress on the intact rock surrounding the fractures, resulting in the closure of them, which in turn increases the P-wave velocity. This phenomenon is in contrast to the classical view, where the propagating fluid pressure diffuses throughout the fracture network and causes a decrease in the P-wave velocity by increasing the pore pressure throughout the fractured medium.

**Abbreviations**

EGS Enhanced geothermal system

DEM Discrete element method

DITE Drilling induced tensile fractures

THM Thermo-hydro-mechanical

WHP Well-head pressure

DHP Downhole pressure

PBOT Pressure build-up onset time

LHS Left-hand side

RHS Right-hand side

Vp P-wave velocity

**Acknowledgements**

We would also like to thank Julie Maury, Thibault Candel, Guido Blöcher and Christophe Fond for their constructive comments. We also thank Hamid Pourpak from Total Energies, Rima Ghazal and Jim Hazzard of Itasca S.A.S. and Itasca C.G. for their technical assistance and fruitful discussions. This work is dedicated to the memory of François Cornet.

**Author contributions**

DJ: methodology, writing—original draft preparation, data curation, and visualization. JS: conceptualization, supervision, and revision. FC: conceptualization and supervision.

**Funding**

This work of the Interdisciplinary Thematic Institute GeoT, 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 SFRI-STRAT'US project (ANR-20-SFRI-001) under the framework of the French Investments for the Future Program, and by Total Energies (Grant #FR00038306).

**Availability of data and materials**

The data sets used and/or analyzed in the current study are available from the corresponding author upon reasonable request.

**Declarations**

**Competing interests**

The authors declare that they have no competing interests.

Received: 5 February 2024 Accepted: 10 December 2024

Published online: 19 February 2025

**References**

Anderson EW. The dynamics of faulting. Trans Edinburgh Geol Soc. 1905;8:387–402.

Ball P. Water: Water—an enduring mystery. Nature. 2008;452:291–2.

Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int J Rock Mech. 1983;20:249–68.

Barth A, Wenzel F, Langenbruch C. Probability of earthquake occurrence and magnitude estimation in the post shut-in phase of geothermal projects. J Seismol. 2013;17:5–11.

Barton N, Bandis S, Bakhtar K. Strength, deformation and conductivity coupling of rock joints. Int J Rock Mech Min Sci Geomech Abs. 1985;22(3):121–40.

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

Calo M, Dorbath C, Cornet FH, Cuenot N. Large-scale aseismic motion identified through 4-D P-wave tomography. Geophys J Int. 2011;186:1295–314. https://doi.org/10.1111/j.1365-246X.2011.05108.x.

Cappa F, Scuderi MM, Collettini C, Guglielmi Y, Avouac JP. Stabilization of fault slip by fluid injection in the laboratory and in situ. Sci Adv. 2019;5:1.

Cauchie L, Lengliné O, Schmittbuhl J. Seismic asperity size evolution during fluid injection: case study of the 1993 Soultz-sous-Forëts injection. Geophys J Int. 2020;221(2):968–80 https://doi.org/10.1093/gji/ggaa051.
Cook NG. Natural joints in rock: Mechanical, hydraulic and seismic behaviour and properties under normal stress. Int J Rock Mech Mining Sci Geomech Abst. 1992;29(3):198-223.

Cornet FH. Elements of Crustal Geomechanics. Cambridge University Press; 2015.

Cornet FH, Bérard Th, Bourouis S. How close to failure is a granite rock mass at a 5km depth? Int J Rock Mech Min Sci. 2007;44:47-66.

Cornet FH. Seismic and aseismic motions generated by fluid injections Geomech. Energy Environ 2016;542-54 https://doi.org/10.1016/j.gete.2015.12.003.

Cundall PA. Formulation of a three-dimensional distinct element model—Part I. A scheme to detect and represent contacts in asystem composed of many polyhedral blocks. Int J Rock Mech Mining Sci Geomech Abst. 1988;25:107-16.

Darot M, Reuchlé T. Acoustic wave velocity and permeability evolution during pressure cycles on a thermally cracked granite. Int J Rock Mech Min Sci. 2000;37:1019-26.

Deng Q, Blöcher G, Cacace M, Schmittbuhl J. Hydraulic difusivity of a partially open rough fracture. Rock Mech Rock Eng. 2021;54:5493-515.

Deng Q, Schmittbuhl J, Cacace M, Blöcher G. Mechanical stiffness and permeability of a reservoir-scale rough fracture during closure. J Geophys Res. 2024;129(9):e2024JB029001.

Dezayes C, Genter A, Valley B. Structure of the low permeable naturally fractured geothermal reservoir at Soulitz. CR Geosci. 2010;342:517-30.

Dublanchet P. Fluid driven shear cracks on a strengthening rate-and-state frictional fault. J Mech Phys Soli. 2019;132:1.

Evans KF. Permeability creation and damage due to massive fluid injections into granite at 3.5 km at Soulitz: 2. Critical stress and fracture strength. J Geophys Res Solid Earth. 2005;110:1.

Fine RA, Millero FJ. Compressibility of water as a function of temperature and pressure. J Chem Phys. 1973;59:5529-36.

Fletcher JB, Sykes LR. Earthquakes related to hydraulic mining and natural seismic activity in western New York state. J Geophys Res. 1977;82:3767-80.

Genter A, Evans K, Cuenot N, Fritsch D, Sanjuan B. Contribution of the exploration of deep crystalline fractured reservoir of Soulitz to the knowledge of enhanced geothermal systems (EGS). CR GEOSCI 2010;342(7-8) https://doi.org/10.1016/j.crte.2010.01.006.

Gentier, S., Genter, A., Sanjuan, B., Bourgine, B., Dezayes, C., Hosni, A., Breton, J.P., Nicol, N., Quinquis, J.P., Crouzet, C., Braibant, G., Brach, M., Moussay, A., Foucher, J.C., Jouin, F., 2003. "Hot Dry Rock" Project; Soulitz-sous-Forêts (No. RP-52730-FR). BRGM, France.

Gentier S, Rachez X, Dezayes C, Blaisonneau A, Genter A. How to understand the effect of the hydraulic stimulation in terms of hydro-mechanical behavior at Soulitz-sous-Forêts (France). GRC Transactions. 2005;29:1.

Haenel, R., 1983. Geothermal investigations in the Rhenish Massif, in: Fuchs, K., von Gehlen, K., Mälzer, H., Murawski, H., Semmel, A. (Eds) Plateau Uplift. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-69219-2_29

Hart RD. An introduction to Distinct Element Modeling for Rock Engineering. Pergamon Press; 1993.

Itasca CG, Inc. Minnesota, USA: Itasca Consulting Group, Inc.; 2016.

Jaeger, J.C., Cook, N.G.W., Zimmerman, R.W., 2007. Fundamentals of rock mechanics, 4th Edition. ed. Blackwell.

Jin L, Zoback MD. Fully coupled nonlinear fluid flow and poroelasticity in arbitrarily fractured porous media: A hybrid-dimensional computational model. J Geophys Res Solid Earth. 2017;122:7626-58.

Kee G, Rummel F. Hydrofrac stress data for the European HDR research project test site Soulitz-sous-Forêts. Int J Rock Mech Min Sci Geomech Abst. 1993;30:973-6.

Kestin J, Sokolov M, Wakeham WA. Viscosity of liquid water in the range -8 °C to 150 °C. J Phys Chem Ref Data. 1978;7:941-8.

Kranz RL, Frankel AD, Engelder T, Scholz CH. The Permeability of Whole and Jointed Barre Granite. Int J Rock Mech Mining Sci Geomech Abst. 1979;16:225-34.

Lengline O, Boubacar M, Schmittbuhl J. Seismicity related to the hydraulic stimulation of GRT1, Rittershoffen. France Geophys J Int. 2017;208:1704-15.

Mavko G, Vanorio T. The influence of pore fluids and frequency on apparent effective stress behavior of seismic velocities. Geophysics. 2010;75:N1-7.

Murphy H, Huang C, Dash Z, Zyvoloski G, White A. Semi-analytical solutions for fluid flow in rock joints with pressure-dependent openings. Water Resour Res. 2004;40:1.

Ohtake M. Seismic activity induced by water injection at Matsushiro. Japan J Phys Earth. 1974;22:163-76.

Orlecka-Sikora B, Lasocki S, Kocot J, Szepieniec T, Grasso JR, Garcia-Aristizabal A, et al. An open data infrastructure for the study of anthropogenic hazards linked to georesource exploitation. Sci Data [Internet]. 2020;7(1):89. https://doi.org/10.1038/s41597-020-0429-3.

Parotidis M, Shapiro SA, Rothert E. Back front of seismicity induced after termination of borehole fluid injection. Geophys Res Lett. 2004;31:1.

Rice JR. Fault stress states, pore pressure distributions, and the weakness of the san andreas fault. Int Geophys. 1992;51:475-503.

Rozhko AY. Role of seepage forces on seismicity triggering. J Geophys Res Solid Earth. 2010;115:1.

Rummel F. Physical properties of the rock in the granitic section of Borehole GPK1. Soulitz-Sous-Forêts Geother Sci Tech. 1991;3:199-216.

Sausse J. Hydromechanical properties and alteration of natural fracture surfaces in the Soulitz granite (Bas-Rhin, France). Tectonophysics. 2002;348:169-85.

Sausse J, Dezayes C, Dorbath L, Genter A, Place J. 3D model of fracture zones at Soulitz-sous-Forêts based on geological data, image logs, induced microseismicity and vertical seismic profiles. CR Geosci. 2010;342:531-45.

Shapiro SA, Huenges E, Borm G. Estimating the crust permeability from fluid-injection-induced seismic emission at the KTB site. Geophys J Int. 1997;131:F15-8.

Shapiro SA, Audigane P, Royer JJ. Large-scale in situ permeability tensor of rocks from induced microseismicity. Geophys J Int. 1999;137:207-13.

Shapiro SA. Fluid-Induced Seismicity. Cambridge University Press; 2015 https://doi.org/10.1017/CBO9781139051132.
Shapiro SA, Rothert E, Rath V, Rindschwentner J. Characterization of fluid transport properties of reservoirs using induced microseismicity. Geophysics. 2002;67:212-20.

Simpson DW, Leith WS, Scholz CH. Two types of reservoir induced seismicity. Bull Seism Soc Am. 1988;78:2025-40.

Stefansson I, Varela J, Keilegavlen E, Berre I. Flexible and rigorous numerical modelling of multiphysics processes in fractured porous media using PorePy. Results Appl Math. 2024; 21100428. https://doi.org/10.1016/j.rinam.2023.100428.

Talwani P, Acree S. Pore pressure diffusion and the mechanism of reservoir-induced seismicity. Pageoph. 1985;122:947-65.

Todd T, Simmons G. Effect of Pore Pressure on the Velocity of Compressional Waves in Low-Porosity Rocks. J Geophys Res. 1972;77:3731-43.

Ucar E, Berre I, Keilegavlen E. Postinjection normal closure of fractures as a mechanism for induced seismicity. Geophys Res Lett. 2017;44:9598-606.

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

Velley B, Evans KF. Stress state at Soultz-sous-Forêts to 5 km depth from wellbore failure and hydraulic observations. Presented at the Workshop on Geothermal Reservoir Engineering: Stanford University, Stanford, California, USA; 2007.

Walsh JB. Effect of Pore Pressure and Confining Pressure on Fracture Permeability. Int J Rock Mech Mining Sci Geomech Abst. 1981;18:429-35.

Wang HF. Theory of linear poroelasticity with applications to geomechanics and hydrogeology. Princeton: Princeton University Press; 2000.

Weidler, R., Gerard, A., Baria, R., Baumgartner, J., Jung, R., 2002. Hydraulic and micro-seismic results of massive stimulation test at 5 km depth at the European Hot-Dry-Rock test site Soultz, France. Presented at the 27th Workshop on Geothermal Reservoir Engineering, Stanford, California, USA, pp. 95-100.

William L., Ellsworth Injection-Induced Earthquakes Movers and Shakers Science 2013;341(6142) https://doi.org/10.1126/science.1225942.

Witherspoon PA, Wang JSY, Iwai K, Gale JE. Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour Res. 1980;16:1016-24.

Yin Z, Huang H, Zhang F, Zhang L, Maxwell S. Three-dimensional distinct element modeling of fault reactivation and induced seismicity due to hydraulic fracturing injection and backflow. J Rock Mech Geothec Eng. 2020;12(752-767):45.

Zangerl C, Evans KF, Eberhardt E, Loew S. Normal stiffness of fractures in granitic rock: a compilation of laboratory and in-situ experiments. Int J Rock Mech Min Sci. 2008;45(8):1500-7.

Zoback MD. Reservoir Geomechanics Cambridge University Press; 2007. https://doi.org/10.1017/CBO9780511586477.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.