# B-Value Evolution During the Strasbourg 2018–2022 Sequence of Induced Earthquakes

Riccardo Minetto, Olivier Lengliné, Jean Schmittbuhl

► To cite this version:

Riccardo Minetto, Olivier Lengliné, Jean Schmittbuhl. B-Value Evolution During the Strasbourg 2018–2022 Sequence of Induced Earthquakes. Journal of Geophysical Research: Solid Earth, 2025, 130 (12), 10.1029/2025JB031744. http://dx.doi.org/10.5194/jes-12458014

HAL Id: hal-05458014

https://univoak.hal.science/hal-05458014v1

Submitted on 14 Jan 2026

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 CC BY 4.0 - Attribution - International License
# JGR Solid Earth

## RESEARCH ARTICLE
10.1029/2025JB031744

**Key Points:**

* Monitoring b-value variations over a 4-year sequence of induced and triggered seismicity

* Clustering analysis of the seismicity reveals spatial clusters with distinct b-values

* Largest magnitude events are associated with low b-value clusters

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

R. Minetto, minetto.riccardo@gmail.com

**Citation:**

Minetto, R., Lengliné, O., & Schmittbuhl, J. (2025). B-value evolution during the Strasbourg 2018–2022 sequence of induced earthquakes. *Journal of Geophysical Research: Solid Earth*, 130, e2025JB031744. https://doi.org/10.1029/2025JB031744

Received 7 APR 2025
Accepted 30 NOV 2025

**Author Contributions:**

**Conceptualization:** R. Minetto, O. Lengliné, J. Schmittbuhl

**Data curation:** J. Schmittbuhl

**Formal analysis:** R. Minetto

**Funding acquisition:** O. Lengliné, J. Schmittbuhl

**Investigation:** R. Minetto

**Methodology:** R. Minetto, O. Lengliné

**Project administration:** O. Lengliné, J. Schmittbuhl

**Software:** R. Minetto, O. Lengliné

**Supervision:** O. Lengliné, J. Schmittbuhl

**Validation:** O. Lengliné, J. Schmittbuhl

**Visualization:** R. Minetto, O. Lengliné

**Writing – original draft:** R. Minetto

**Writing – review & editing:** O. Lengliné, J. Schmittbuhl

© 2025. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# B-Value Evolution During the Strasbourg 2018–2022 Sequence of Induced Earthquakes

R. Minetto¹ ID, O. Lengliné¹ ID, and J. Schmittbuhl¹ ID

¹EOST/ITES, Strasbourg University/CNRS, Strasbourg, France

**Abstract** The b-value, defined as the slope of the cumulative earthquake frequency–magnitude distribution in semi-logarithmic space, quantifies the relative occurrence of small and large earthquakes. It is a key metric for hazard assessment and is often used in adaptive traffic light systems to mitigate induced seismicity. In this study, we analyze a 4-year sequence of induced earthquakes at the Vendenheim geothermal site near Strasbourg (France), where multiple sequential operations were carried out near a regional fault zone, to explore the spatio-temporal variability of the b-value during fault activity. Our analysis highlights the importance of spatial clustering in understanding b-value variability. Despite occurring in close spatial and temporal proximity, seismicity clusters exhibit markedly different b-values. The clusters with the lowest b-values host the largest magnitude events ($3 < M_L < 4$). Even when significant seismic migration is observed within the clusters, their b-value remains relatively constant over time. The observed temporal fluctuations in the b-value of the entire induced seismicity population are shown to be closely related to the relative activity of each cluster over time. We propose that b-value variations, along with the observed seismicity migration, are related to inherited spatial variability in the geomechanical properties of faults rather than to the dynamic evolution of stress conditions.

**Plain Language Summary** The relative frequency of small to large earthquakes in a given region, known as the b-value, is a key parameter for estimating earthquake recurrence times and assessing the probability of large earthquakes. Its estimation is particularly relevant in contexts where human activities, such as fluid injections for geothermal energy or hydrocarbon extraction, have the potential to induce large earthquakes. However, the b-value is generally found to vary across space and time, and identifying the factors controlling these variations remains challenging. This study examines a 4-year period of induced seismicity associated with geothermal energy operations in Vendenheim, France. We demonstrate that spatial and temporal variations in the b-value result from the activation of multiple fault regions, each exhibiting distinct b-values. We attribute these differences to variations in conditions along the faults, such as material properties and pre-existing stress. Our findings highlight the importance of accurately identifying the activated geological structures to better assess the seismic hazard associated with human activities.

## 1. Introduction

The magnitude distribution of earthquakes has long been a classical tool for characterizing seismicity. This distribution typically follows an empirical relationship known as the Gutenberg-Richter law:

$$ \log_{10}N(M) = a - bM \quad (1) $$

where $N(M)$ is the cumulative number of earthquakes with magnitude larger than or equal to M (Gutenberg & Richter, 1944). The b-value provides information on the ratio between small and large events, while a represents the total number of earthquakes with $M \ge 0$, serving as a measure of earthquake productivity.

In the context of seismic sequences, the b-value is generally analyzed either over time (e.g., Herrmann et al., 2019) or in space within regions of high event density (e.g., Igonin et al., 2018). The interpretation of the spatiotemporal evolution of the b-value has received considerable interest over the past decades (El-Isa & Eaton, 2014; Herrmann et al., 2022). Numerous studies have shown that local variations in the b-value reflect an inverse relationship with the stress levels of the activated region. Differential stress is most commonly cited as the primary type of stress controlling the b-value (Scholz, 2015). Several other relationships, such as those between b-value and depth (Spada et al., 2013), Coulomb stress (Gulia et al., 2018; Scholz, 1968), and stress-drop (Wyss
et al., 1973), are thought to reflect this dependency (Scholz, 2015). However, studies suggest that other types of stress can directly control the b-value, such as shear (Mukuhira et al., 2024) and normal stress (Dublanchet, 2020).

In addition to stress conditions, some studies suggest that the b-value is influenced by fault properties, such as the size of the fault system (e.g., Chen et al., 2018; Kettlety et al., 2020), fault roughness (Goebel et al., 2017), strength heterogeneities along the fault (Steacy et al., 1996), lithology (Collettini & Tinti, 2025), and stiffness (Stormo et al., 2015).

An relevant context for studying the b-value evolution is the injection of fluids into the subsurface. This technique is widely used to improve the permeability of rocks, enabling applications such as the extraction of shale gas, wastewater disposal, or the production of geothermal energy (Grigoli et al., 2017; Keranen & Weingarten, 2018; Norris et al., 2016). This process alters the stress state of the subsurface through mechanisms such as pore pressure diffusion (Shapiro et al., 1997), aseismic slip (De Barros et al., 2020), and poroelastic effects (Ge & Saar, 2022). Past injection operations have shown that such stress changes can induce large earthquakes, leading to public concern, infrastructure damage, and the suspension of ongoing and future operations (Candela et al., 2018; Keranen & Weingarten, 2018, and references therein). Identifying regions with a higher probability of hosting large earthquakes during planned injections is therefore critical to mitigate these risks (Muntendam-Bos & Grobbe, 2022).

The Geoven deep geothermal energy site at Vendenheim, France, offers a valuable opportunity to explore b-value variations during fluid injection due to its long operational period (2018–2021), multiple injection phases, and a geological setting with numerous preexisting, critically stressed faults. A first study by Schmittbuhl et al. (2021) provided a high-resolution earthquake catalog of the induced seismic sequence up to the shut-in of the wells in January 2021. Their analysis focused on exploring the relationship between seismicity, injection operations, and known geological structures. A second study by Lengliné et al. (2023) introduced a 3D velocity model of the area and analyzed the source properties of the four main events, including the largest event ($M_L$3.9) that occurred on 26 June 2021, six months after shut-in.

Here, we use spatial clustering to explore key aspects relevant to hazard analysis, such as the imaging of newly active faults, the space-time variability of the b-value, and seismicity migration. To achieve this, we created a new high-resolution catalog for the entire induced seismic sequence associated with the development of the deep geothermal reservoir, up to October 2023. Using this catalog, we show that cluster analysis is essential for uncovering spatial variations in the b-value and propose that these variations, along with the observed seismicity migration, are related to inherited spatial variability in the geomechanical properties of faults rather than differences or time evolution of stress conditions.

The article is structured as follows: in Section 2, we introduce the waveform database, describe the creation of a new catalog based on template matching and double-difference relocation, and explain the approach for seismic cluster identification. In Section 3, we present the cluster features in terms of geometry, Coulomb stress, frequency-magnitude distribution, and space-time migration. In Section 4, we discuss the cluster features in relation to the injection operation and fault properties, and propose potential causes for the variations in the b-value. In Section 5, we formulate the main conclusions of our study.

## 2. Material and Methods

### 2.1. Waveform Database

The seismicity catalog provided by the French seismological survey BCSF-RéNaSS (Bureau Central Sismologique Français and Réseau National de Surveillance Sismique) includes 582 events that occurred between January 2018 and September 2023 within 10 km of the injection site. These events are distributed across two main clusters: a northern cluster in the vicinity of the wells and a southern cluster approximately 5 km south of the wells. The BCSF-RéNaSS catalog was created using SeisComP (GFZ & GEMPA GmbH, 2008) and the contribution of human operators who visually detected earthquakes, inspected automatic detections, and refined wave arrival picks.

For this study, we considered only the stations that contributed the most to seismicity detection by selecting those with at least 50 P-wave picks in the BCSF-RéNaSS catalog. This is a total of 38 stations (Figure 1 and Figure S1 in Supporting Information S1). Of these, 11 are Raspberry Shake® stations, short-period stations with one or three
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Map of the study area. Red inverted triangles indicate Raspberry Shake® stations, while blue inverted triangles represent broadband stations. Black circles mark the relocated epicenters of the earthquakes that occurred between January 2018 and October 2023. Seismicity is overall structured in two groups of clusters: northern clusters near the wells and southern clusters around the $M_L$3.0 event. Green lines indicate the simplified trajectories of the GT1 and GT2 wells. Light brown areas represent the projections of known fault zones between 3 and 6 km depth, as derived from the regional GEORG geological model. The focal mechanism of the four largest events are also shown. See Figure S1 in Supporting Information S1 for the locations of all stations used in this study.</figcaption></figure>

components recording at 50 and 100 Hz, respectively. These stations have the unique feature of being hosted by citizens in their homes. The remaining 27 stations are broadband stations or accelerometers with a sampling frequency of either 100 or 200 Hz. 13 of these stations belong to the French permanent network, 6 to Fonroche Géothermie, 3 to a local network installed during the Vendenheim induced sequence, 2 to the French accelerometric network, 2 to the Erdbebendienst Südwest network, and 1 to the German regional seismic network. We point out that some of the selected stations were not active for the entire period under investigation (January 2018–October 2023). As a result, the number of active stations at any given time can be highly variable (Figure S2 in Supporting Information S1).

## 2.2. Catalog Creation

### 2.2.1. Template Matching

Template matching, or matched filtering, is a technique used to maximize the number of detected earthquakes (Gibbons & Ringdal, 2006). It consists of cross-correlating a set of reference events, called templates, with continuous recordings to extract earthquakes with similar waveforms. Template matching has been widely used in different geological settings to study natural and induced earthquakes (e.g., Cochran et al., 2018; Shelly et al., 2016).

We used the events of the BCSF-RéNaSS catalog as templates, keeping only those with a signal-to-noise ratio (SNR), averaged over the three components, greater than three at a minimum of three stations (458 out of 582 events). This step helps reduce the computational cost of template matching and removes low-quality waveforms that could cause false detections. Each template is 5.12 s long and starts 0.3 s before the P-wave arrival. Both templates and continuous recordings were resampled at 100 Hz and filtered between 8 and 30 Hz.

Template matching was performed by cross-correlating the seven highest-quality components of each template with the corresponding 24-hr long traces. Seven components were chosen because this number maximizes the number of detections while maintaining a low rate of false positives. The seven components were selected as follows. First, the components available for the current day were selected. The components of each station were then grouped together, with the vertical component (Z) ranked first, followed by the horizontal components (N and E). Finally, the different groups of components were sorted according to their average SNR and the first seven components were selected.

A detection is issued when the stack of the seven cross-correlation functions (time-shifted to take into account the travel time differences between stations) exceeds a template-based correlation coefficient (CC) threshold. This threshold is estimated by performing the cross-correlation using a flipped and reversed version of the seven selected components. From the resulting stacked correlation function, we derived the cumulative distribution of the CC values. The CC threshold is then defined as the CC for which the cumulative distribution is above a given false alarm rate. The advantage of a template-based CC threshold is that it takes into account the temporal variation of noise conditions and, more generally, the variability of the information contained in each template. The time tolerance between two consecutive events was set to 4.5 s.

Certain stations can cause false detections due to sensor glitches, which produce waveforms that are highly similar to each other and resemble impulsive P-wave arrivals. As a result, these glitches can raise the CC above the detection threshold. To address this issue, we first verified whether the station producing the highest average CC with the template is among those known to experience glitches. If so, we cross-correlated its components with representative glitch waveforms from that station. If the CC for the glitch was higher than that for the template, we discarded the detection.
We detected 2,333 events after template matching, including the templates. This is a four-fold increase compared to the starting catalog. As a result, the global magnitude of completeness decreases from 0.75 to 0.4 (Figure S3a in Supporting Information S1).

### 2.2.2. Magnitude Estimation

We estimated the magnitude of each detection based on the magnitude of the template that detected it. For each component with a CC greater than 0.35, we first applied a band-pass filter between 8 and 24 Hz. This frequency band was selected to ensure that the signal remained below the corner frequency (as 99% of the detections are associated with templates with magnitudes below 1.8) while simultaneously maximizing the SNR of the smallest events. We then calculated the local magnitude difference ($\Delta M_L$) between the detection and the template using the following relation:

$$ \Delta M_L = \frac{\log_{10} \frac{A_T}{A_D}}{1.143} \qquad (2) $$

where $A_D$ and $A_T$ are the peak-to-peak amplitudes of the detection and the template, respectively, and 1.143 is the coefficient of the equation $M_0 = 10^{1.143M_L+9.81}$, obtained by linear regression of the seismic moment ($M_0$) versus $M_L$ for the largest events (Figure S4 in Supporting Information S1). We then calculated a weighted average of the $\Delta M_L$ values, using the CC values as weights to emphasize the components with higher similarity. Finally, the resulting value was added to the magnitude of the template to derive the magnitude of the detection.

### 2.2.3. Event Relocation

The earthquakes of the BCSF-RéNaSS catalog were located using a 1D velocity model based on Rothé and Peterschmitt (1950). To obtain more accurate absolute locations, we located the templates using the local 3D velocity model proposed by Lengliné et al. (2023) and the NonLinLoc software (Lomax et al., 2000). Each detection was then initially assigned the location of its corresponding template, and P- and S-wave arrival times for the detection were derived by assuming identical travel times to those of the template. With this new information, we estimated the differential times between pairs of earthquakes observed at common stations. These values were obtained by cross-correlating waveforms filtered between 8 and 24 Hz, lasting 0.8 s for P waves and 1 s for S waves, and starting 0.2 s before the wave arrival. Only differential times with CC > 0.65 and pairs with at least 9 differential times were kept to ensure reliable measurements.

The differential times were then used to relocate the events with a double-difference method (Trugman & Shearer, 2017). Double-difference methods are able to provide high-resolution relative locations by taking advantage of the fact that closely located events share a similar path, which helps to minimize errors due to unmodeled velocity structures. During the relocation process, we assigned a greater weight to differential times calculated from template events (i.e., the events with directly measured absolute locations) to better constrain the relocation process. Finally, the location uncertainties were estimated using a bootstrap approach with 2,000 resamples (Trugman & Shearer, 2017).

Of the 2,333 detections, 1,270 could be relocated. The average location errors calculated with bootstrap are 80 and 180 m for the horizontal and vertical directions, respectively. We point out that location uncertainty is time-dependent, as the number of active stations in the area changes over time (Figure S2 in Supporting Information S1). For example, from the end of 2019 to the end of the study period, the average horizontal and vertical uncertainties decrease to 50 and 120 m, respectively. Figure S5 in Supporting Information S1 provides a comparison of these new locations with those from the BCSF-RéNaSS catalog and the catalog of Schmittbuhl et al. (2021).

### 2.3. Cluster Identification

We define a cluster as a group of closely spaced events that align along a common plane in 3D space. Clusters with different orientations or varying earthquake density can occur in close proximity to one another. This proximity makes it difficult to separate clusters using a single distance threshold, as required by algorithms such as DBScan (Density-Based Spatial clustering) (Ester et al., 1996). To address this limitation, we combined DBScan with an expectation-maximization (EM) algorithm, which fits the data to a specified number of Gaussian
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Map view of the northern and southern clusters and their associated FMD. (a) Map view of the northern clusters. Earthquakes are shown as circles scaled by magnitude. Each cluster is assigned a unique color, while unclustered earthquakes are represented by empty circles. (b) Same as panel a, but for the Southern clusters. (c) Cumulative FMDs of the northern clusters. The black lines represent the GR law estimated with the method of Aki (1965). (d) Like panel c, but for the southern clusters. Each cluster is labeled with a letter ("N" for northern clusters and "S" for southern clusters) followed by a number indicating its order of occurrence. For a depth view of the earthquakes, see Figures S7 and S8 in Supporting Information S1. For a comparison between the b-value calculated using the methods of Ogata and Katsura (1993) and of van der Elst (2021), see Figure S11 in Supporting Information S1.</figcaption></figure>

distributions (Reynolds, 2009). The advantage of this algorithm is that it is sensitive to the planar alignment of earthquakes.

To perform the clustering, we first divided the study area into four subregions, corresponding to the regions with the highest earthquake density (Figure S6 in Supporting Information S1). For regions 1 and 2, we could extract the major clusters by simply applying DBScan with a distance threshold equal to the average location uncertainty of the events (Figures S6a and S6b in Supporting Information S1). Clusters with fewer than 50 events or without clear planar alignment were removed. For region 3, we initially applied DBScan, as we did for regions 1 and 2, to remove isolated events (Figure S6c in Supporting Information S1). We then used EM, setting the number of clusters to two, which corresponds to the number of clusters that we visually identified (Figure S6d in Supporting Information S1).

Region 4 is the most complex due to a NNE-SSW-oriented structure overlapping multiple NW-SE-oriented structures (Figure S6e in Supporting Information S1). We first applied DBScan as in the previous regions. Next, we used EM, gradually increasing the number of clusters until the NNE-SSW structure was clearly isolated, which required a total of five clusters (Figure S6d in Supporting Information S1). Once this structure was isolated, we re-applied DBScan to cluster the remaining events associated with the NW-SE-oriented structures (Figure S6f in Supporting Information S1). We point out that cluster identification is performed using the relocated earthquakes. Earthquakes that could not be relocated were assigned to the same cluster as the template event that detected them.

With this process, we identified eight main clusters: four in the north (N1, N2, N3, N4) and four in the south (S1, S2, S3, S4) of the study region (Figures 2a and 2b, Figure S7 and S8 in Supporting Information S1). Together,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Maps showing the epicenters of the events occurred during (red) and after (blue) each injection phase. (a) 28 February 2018–27 April 2019. (b) 27 April 2019–25 June 2019. (c) 25 June 2019–29 July 2019. (d) 29 July 2019–3 October 2019. (e) 3 October 2019–5 November 2019. (f) 5 November 2019–19 August 2020. (g) 19 August 2020–16 September 2020. (h) 16 September 2020–5 November 2022. Post-injection seismicity includes events occurring from the end of the injection to the start of the next injection. For reference, the events that do not belong to the analyzed period are represented by white circles. The black lines represent the trajectories of the wells GT1 and GT2 from which fluids are injected, while the gray lines the trajectories of the inactive wells.</figcaption></figure>

these clusters account for about 60% of the detected events and include all the largest events of the sequence ($M_L > 2$) (Figure S3b in Supporting Information S1). All of them, with the exception of cluster N2, are located within the granitic basement (Figure S9 in Supporting Information S1). We note that the high-resolution locations obtained through the relocation step were essential for identifying these clusters.

## 2.4. Time Evolution

Injection operations at the Geoven site took place between 2018 and 2021 along two wells (GT1 and GT2, see Figure 3), both reaching a maximum depth of about 5 km (Terrier et al., 2022). The wells targeted the Robertsau Fault (Figure 1), specifically at the transition between the Hercynian granitic basement and the sedimentary cover (Schmittbuhl et al., 2021; Terrier et al., 2022). Fluids were injected from GT1 and GT2 during four and five phases, respectively (Figure 3).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Stress conditions along the fault planes determined from the seismicity clusters: (a) Polar plot showing the variation of the Coulomb stress with fault dip and strike, supposing a friction coefficient of 0.6, a cohesion of 1 MPa, and a depth of 4,300 m (mean depth of all the clustered earthquakes). Green diamonds mark the fault plane orientations determined from the clusters, red diamonds indicate fault planes orientations derived from focal mechanisms, and light brown diamonds the orientation of the faults found in the GEORG model and shown in Figure 1. Strike follows the right-hand rule (dip direction = strike + 90°). (b) Coulomb stress for each cluster, calculated for two different friction coefficients (0.6 and 0.8). (c) Normal, shear, and differential stress at each cluster location. The values in panels (b, c) were estimated supposing a fault depth equal to the mean depth of the corresponding cluster. Shaded bands represent the 1σ uncertainty estimated from the Monte Carlo approach.</figcaption></figure>

<table><caption>**Table 1**

Summary of the Geometry, Duration, and Size of Each Cluster</caption><thead><tr><th>Cluster ID</th><th>Strike (°)</th><th>Dip (°)</th><th>Barycente (lat., lon., depth)</th><th>Extension (along strike)</th><th>Extension (depth)</th><th>Start time</th><th>End time</th><th>N. events (detected, relocated)</th></tr></thead><tbody><tr><td>N1</td><td>N25°E [13°, 37°]</td><td>83°NW [77°, 90°]</td><td>48.6507, 7.7972, 3.9</td><td>0.3</td><td>0.4</td><td>2018-03-14</td><td>2020-08-30</td><td>219, 94</td></tr><tr><td>N2</td><td>N27°E [20°, 35°]</td><td>87°SE [83°, 90°]</td><td>48.6578, 7.7938, 3.3</td><td>0.4</td><td>0.7</td><td>2019-05-06</td><td>2020-12-06</td><td>274, 128</td></tr><tr><td>N3</td><td>N148°E [137°, 158°]</td><td>89°NE [82°, 90°]</td><td>48.6575, 7.8057, 4.6</td><td>0.3</td><td>0.4</td><td>2019-08-28</td><td>2020-12-12</td><td>140, 61</td></tr><tr><td>N4</td><td>N130°E [125°, 135°]</td><td>84°NE [80°, 87°]</td><td>48.6568, 7.8098, 5.0</td><td>0.7</td><td>1.1</td><td>2020-02-22</td><td>2022-02-07</td><td>208, 148</td></tr><tr><td>S1</td><td>N13°E [6°, 20°]</td><td>88°SE [82°, 90°]</td><td>48.6116, 7.7926, 3.9</td><td>0.4</td><td>0.6</td><td>2019-11-11</td><td>2019-12-08</td><td>94, 54</td></tr><tr><td>S2</td><td>N119°E [102°, 135°]</td><td>80°NE [74°, 90°]</td><td>48.6090, 7.7931, 4.1</td><td>0.2</td><td>0.7</td><td>2019-11-12</td><td>2020-01-10</td><td>79, 47</td></tr><tr><td>S3</td><td>N140°E [129°, 150°]</td><td>88°NE [84°, 90°]</td><td>48.6143, 7.7917, 4.5</td><td>0.3</td><td>1.0</td><td>2019-11-15</td><td>2021-08-17</td><td>87, 67</td></tr><tr><td>S4</td><td>N126°E [119°, 133°]</td><td>87°NE [80°, 89°]</td><td>48.6118, 7.7908, 4.9</td><td>0.7</td><td>1.3</td><td>2019-11-21</td><td>2021-05-04</td><td>265, 168</td></tr></tbody></table>

Note. The table lists the strike and dip of the best-fitting plane with associated 95% confidence intervals, the cluster barycenter (latitude, longitude, depth in km), the spatial extent along strike and in depth (km), the start and end dates of activity, and the total number of detected and relocated events.
## 3. Cluster Features

### 3.1. Clusters and Fault Plane Orientations

The clusters align along high-angle planes (80–89°) and along two sets of orientations: NNE-SSW and NW-SE (Figure 4a; Table 1; Figure S10 in Supporting Information S1). These orientations are observed in both the northern and southern regions (Figures 2a and 2b).

The orientations obtained from earthquakes with focal mechanism (Schmittbuhl et al., 2021) match the orientation of the cluster to which each event belongs (Figure 4a). We also note that according to the latest 3D geologic model of the French Rhine Graben (Capar, 2024) NNE-SSW-oriented faults are present in the area, whereas NW-SE-oriented faults are absent (Figures 1 and 4a).

### 3.2. Coulomb Stress Estimation

We estimated the Coulomb stress (CS) for each plane delineated by the seismicity clusters. CS is a measure that quantifies the proximity of a fault plane to rupture (e.g., Beeler et al., 2000). It is defined as:

$$CS = \tau - \mu(\sigma_n - P_f) - C_0$$

where $\tau$ is the shear stress, $\sigma_n$ the normal stress, $P_f$ the pore pressure, $\mu$ the friction coefficient, and $C_0$ the cohesion. Positive CS values indicate unstable (or critical) fault conditions (i.e., slip is favored), while negative CS values indicate stable conditions.

To estimate the $\tau$ and $\sigma_n$ acting on each fault plane, we first determined the orientation (strike and dip angle) of the seismicity clusters using Principal Component Analysis (Jolliffe, 2002). To determine the uncertainty in strike and dip of each cluster, we repeated the estimation of the cluster orientations 2,000 times using the bootstrap samples derived from the relocation step and then derived the 95% confidence interval.

We then assumed a uniform stress field equal to the one measured in Soultz-sous-Forêts, the nearest site with available stress measurements (about 40 km from the Geoven site). The principal stresses (in MPa) below a depth of 3,600 m, assuming a constant density of 2,550 kg/m³, are:

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

$$S_h = 0.537S_v$$

$$S_H = 1.08S_v$$

where $z$ is the depth in meters, $S_v$ is the vertical stress, $S_h$ is the minimum horizontal stress, and $S_H$ is the maximum horizontal stress (Evans et al., 2009). The stress regime is strike-slip, with $S_H$ oriented at an azimuth of 169° (Cornet et al., 2007). We then converted the principal-stress tensor into the geographic coordinate system using a change-of-basis (rotation) matrix, and projected the transformed stress tensor onto each fault plane to obtain the normal and shear stress components (Zoback, 2007).

To account for uncertainties in fault geometry and principal stresses when estimating $\tau$, $\sigma_n$, and CS, we applied a Monte Carlo simulation with 2,000 realizations. In each realization, stresses were recomputed after perturbing the input parameters according to their respective uncertainties. The $S_H$ azimuth was perturbed using a wrapped normal distribution (to account for the 0°–360° domain) with a standard deviation of 14°. The vertical stress $S_v$ was perturbed using a normal distribution with a standard deviation of 0.0098z + 0.6 MPa, following Evans et al. (2009). For each cluster, the strike and dip faults were sampled without replacement from the 2,000 strike-dip pairs derived from the 2,000 bootstrap catalogs generated during the relocation procedure.

The confidence intervals for cohesion, pore-pressure gradient, and friction coefficient are either unavailable or poorly constrained at Vendenheim and nearby geothermal sites (Cornet et al., 2007). We therefore estimated CS assuming a fault cohesion of 1 MPa, as proposed by Schmittbuhl et al. (2021), hydrostatic pore pressure (as typically expected in formations without evidence of fluid overpressure), and two friction coefficients (0.6 and 0.8), which span the typical range for rocks of the upper crust (Byerlee, 1978). From the 2,000 Monte Carlo realizations performed for each friction coefficient, we calculated the standard deviation of $\tau$, $\sigma_n$, and CS to
quantify the uncertainty of these stresses. We note that the hydrostatic pore pressure is primarily used as a baseline to facilitate comparison between clusters. In practice, the pore pressure at the faults is expected to increase during injection, with the magnitude of the increase depending on factors such as the distance from the injection point, the duration and rate of injection, and the hydromechanical properties of the reservoir (Sáez & Lecampion, 2023).

Under the assumed conditions, CS analysis indicates that faults matching the cluster orientations are, in most cases, unstable or close to instability (i.e., CS is close to or greater than zero (Figures 4a and 4b). This critical state persists even when higher cohesion values are used, up to approximately 10 MPa under the most favorable conditions (μ = 0.6). Therefore, even small reductions in effective stress due to pore pressure increases from injection may be sufficient to induce fault slip. Regardless of absolute CS magnitudes, the estimated fault orientations exhibit comparable Coulomb stresses and are among the most critically stressed within the regional stress field (Figure 4a).

### 3.3. Frequency-Magnitude Distribution

We estimated the b-value of each cluster using all detected events and three different methods: Aki (1965), Ogata and Katsura (1993), and van der Elst (2021). The method of Aki (1965) is a maximum likelihood estimation of the b-value, which assumes the validity of the Gutenberg-Richter law (Equation 1). When applying this method, we estimated the magnitude of completeness ($M_c$) with the maximum curvature method (Wyss et al., 1999), adding 0.2 as suggested by Woessner and Wiemer (2005). The method of Ogata and Katsura (1993) allows the simultaneous estimation of $M_c$ and the b-value by maximizing a log-likelihood function that takes into account the entire Frequency-Magnitude distribution (FMD). The method of van der Elst (2021) is based instead on magnitude differences and does not require prior knowledge of $M_c$. Finally, we assessed the uncertainties affecting the b-values with bootstrap resampling (Efron & Tibshirani, 1994).

Clusters N4 and S1, which include the largest magnitude events, have a b-value of approximately 0.6. This is significantly lower than the b-value of the other clusters, which is close to or greater than 1 (Figure 2 and Figure S11 in Supporting Information S1). We note that the b-value in clusters N4 and S1 remains low even when removing the largest events ($M_L$ > 3; Figure S12 in Supporting Information S1) and is stable over time (Figure S13 in Supporting Information S1). In addition, the differences in the b-value between clusters are preserved even without the applied magnitude correction (Equation 2) (Figure S14 in Supporting Information S1).

To evaluate the effect of clustering on the b-value, we calculated the b-value over larger areas, corresponding to the four main regions from which the main clusters were originally derived (Figure S6 in Supporting Information S1). We also computed the temporal evolution of the b-value during and after the last injection period (October 2020-May 2021), when three clusters occurred simultaneously. The temporal analysis is based on the method of van der Elst (2021), which, among the three methods we used, is the least affected by short-term variations in $M_c$. We applied this method using a sliding window of 110 events, shifted by 10 events at each step. This window size corresponds to the smallest window needed to consistently have more than 50 positive magnitude differences to estimate the b-value at each iteration.

If the b-value is calculated over time without any clustering, low b-values appear much later in the sequence, when seismicity is mainly composed of events from the low b-value cluster (Figure 5). Similarly, a crude clustering (Figure S15 in Supporting Information S1) completely obscures the presence of low b-values in the southern region (Figure S15b in Supporting Information S1). In the northern region, a low b-value may still be identified, but the FMD is highly bilinear, making the interpretation of the b-value ambiguous (Figure S15d in Supporting Information S1).

### 3.4. Relation Between Clusters and Injection Operations

Cluster N1 is the first of the main clusters to be activated, and its occurrence is highly correlated with injections along the GT1 well (Figure 6 and Figure S7 in Supporting Information S1). However, this cluster accounts for only about 25% of the seismicity induced by the GT1 injections, which mainly consists of small clusters of varying shapes and orientations (Figure S7 in Supporting Information S1).

The occurrence of clusters N2, N3, and N4 is closely associated with injections along GT2 (Figure 6 and Figure S8 in Supporting Information S1), and 80% of all seismicity occurring during GT2 injections belongs to these clusters. N2 activated during three of the five GT2 injection phases, N3 during two, while cluster N4—the one that
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Temporal evolution of the b-value and its association with the activated clusters (a) Variation of the b-value with time during and immediately after the last injection. The shaded area represents the 95% confidence interval. The time of the measurement corresponds to the average time of the events in the moving window. (b) Percentage of the events belonging to each of the identified clusters active in the chosen time period. For reference, the b-value of each active cluster is reported in the legend.</figcaption></figure>

hosted the three largest events of the sequence (3.3 < *M*<sub>L</sub> < 3.9)—activated during the last injection (Figure S7 in Supporting Information S1). In particular, cluster N4 differs from the other northern clusters in that it remained active for more than a year after the end of injection (shut-in).

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Temporal evolution of the relocated seismicity and relation with the injection operations. (a) Latitude variation of earthquakes over time. (b) Hypocentral depth variation over time. The gray and light gray shaded areas indicate periods when well operations (hydraulic tests or circulation tests) were performed at the GT1 and GT2 wells, respectively. Each circle represents an event, with size scaled by magnitude and color corresponding to the cluster it belongs to.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

**Figure 7.** Relation between bursts (i.e., periods of increased seismicity rate) and evolution of the seismicity within clusters with different b-value. (a) Temporal evolution of the normalized seismicity rate (black bars) of a low b-value cluster (N4); bursts are highlighted by shaded areas of different colors. All detected events with magnitude greater than the magnitude of completeness of the cluster are considered. (b) Same as panel a, but for a high b-value cluster (S4) (c) along-strike cross section, with each polygon outlining the area occupied by the events (points) occurring during a burst of cluster N4 and shown in panel (a) Polygons and earthquakes associated to the same burst follow the same color coding as in panel (a) Only relocated events having an average distance from other cluster events within one standard deviation are used. (d) Same as panel c, but for cluster S4. (e) Burst area (calculated from the polygons shown in panel c and $M_{max}$ of each burst of cluster N4. Panel (f) Same as panel (e), but for cluster S4. Figure S16 in Supporting Information S1 provides the same representation for additional clusters.

The connection between the southern clusters (S1, S2, S3, and S4) and injections is less evident, as the seismicity occurred approximately 4 km south of the closest injection point (Figure 3). Schmittbuhl et al. (2021) showed that the delay between the GT1 injection and the onset of seismicity in the south is about 10 days, which is consistent with fluid pressure migration. This observation was deduced from the triggering of seismicity in the Southern cluster after two GT1 injections with similar durations and flow rates: one in November 2019 and one in late August 2020 (Schmittbuhl et al., 2021).

### 3.5. Migration

All identified clusters show evidence of migration, either downward (clusters S2, S3), upward (cluster N4), or in both downward and upward directions (cluster S4) (Figure 6, Figures S7 and S8 in Supporting Information S1). Migration is observed along both strike and dip, although the latter is the dominant direction (Figures S7 and S9 in Supporting Information S1).

The observed migrations are not a continuous process; instead, they occur intermittently through successive bursts of seismicity (Figures 7a and 7d, and Figure S16 in Supporting Information S1), where a burst is defined as a period of increased seismic activity within the cluster. Although in some cases bursts coincide with the start of injection (clusters N1–N3), in other cases bursts occur without an active injection (S1–S4, and N4) (Figure 6). Bursts tend to predominantly activate new fault regions, leaving previously active regions inactive (Figure S16 in Supporting Information S1). This is particularly evident in clusters that persist after injection (Figures 7c and 7d).

The area activated during a burst varies between clusters. In low b-value clusters (e.g., N4 and S1), the burst area tends to increase with each burst and is larger than in high b-value clusters (e.g., S2–S4), where it remains small over time (Figures 7e and 7f, and Figure S16 in Supporting Information S1). Similarly, the evolution of the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Migration of seismicity within clusters N4 and S4. (a) Variation of distance from the centroid of the first five earthquakes over time for cluster N4. (b) Same as panel a, but for cluster S4. (c) Variation of distance from the centroid of the first five earthquakes with event index for cluster N4. (d) Same as panel c, but for cluster S4. Each point represents an event, colored by its origin time. Shaded bands indicate periods of fluid injection. Migration for the other clusters is shown in Figure S17 of Supporting Information S1.</figcaption></figure>

maximum magnitude ($M_{max}$) after each burst differs between clusters. In low b-value clusters, $M_{max}$ gradually increases with each burst (Figure 7e and Figure S16 in Supporting Information S1) and reaches larger values than in high b-value clusters, where it remains small throughout their entire duration (Figure 7f and Figure S16 in Supporting Information S1).

In addition, clusters that remain active after injection and propagate along a dominant direction (e.g., S2, S3, S4) tend to migrate linearly, although discontinuously, over time (Figure 8 and Figure S17 in Supporting Information S1). When migration is analyzed using the event index, the trend appears linear and continuous in most post-injection clusters (S1–S4, and N4) (Figures 8c and 8d, and Figure S17 in Supporting Information S1).

The migration process can last for the entire duration of the cluster, which can range from a few weeks to about one year (Figure 8 and Figure S17 in Supporting Information S1). By calculating the distance of the cluster events from the centroid of the first five earthquakes and fitting a line to the distances within the 70th to 90th percentiles over a moving window of 15 events, we obtained migration velocities on the order of tens of meters per day for clusters S2–S4 and a few meters per day for cluster N4 (Figure S17 in Supporting Information S1).

In clusters that occurred almost exclusively during injection (N1, N2, N3), we were unable to identify clear general migration patterns due to the complexity of the migration (Figures S17a–S17c and S17e–S17g in Supporting Information S1). Similarly, the location accuracy and the small number of events during the burst do not allow for the determination of clear migration patterns during bursts.

# 4. Discussion

## 4.1. B-Value as a Fault Property

We have shown that seismicity at the Geoven deep geothermal site predominantly clusters along eight fault planes, which we will hereafter refer to as fault segments, as they may be part of a larger fault. Two of the eight fault segments are characterized by significantly lower b-values compared to the others (Figures 2c and 2d). We
have also shown that low b-values can only be detected if seismicity clustering is precisely identified. Indeed, if a temporal analysis is performed during periods when several clusters are active (Figure 5) or if a less precise distance-based clustering is applied (Figure S15 in Supporting Information S1), low b-values either disappear completely or lose their interpretative value. This is because the FMDs can become highly bilinear (i.e., have an excess of large events). Bilinear FMDs, which are typically observed in the region where the largest event occurred (Igonin et al., 2018; Kettlety et al., 2020), may therefore indicate a mixing of distinct FMDs from different fault structures. Our observations thus suggest a strong dependence of the b-value on the fluid-stimulated space and, more specifically, a dependence of the b-value on the local conditions along the fault segments.

## 4.2. Variations in Fault Rupture Behavior

With the exception of cluster S1, the observed fault segments did not activate all at once but in successive bursts of activity (Figures 7a–7d and Figure S16 in Supporting Information S1). Although in some clusters increased seismicity is linked to injection and likely results from stress changes caused by high-pressure fluids, bursts can also occur months after injection. In such cases, seismicity is characterized by a generally linear migration over distance that is discontinuous in time but continuous in event index (Figure 8 and Figure S17 in Supporting Information S1). These migration patterns are consistent with earthquake interactions, a self-sustaining process in which ruptures trigger the nucleation of subsequent ruptures (Fischer & Hainzl, 2021). In this scenario, injected fluids still play a significant role, as they are needed to initiate seismic activity.

We have also shown that $M_{max}$ and burst area remain small over time in fault segments with high b-values, whereas they become large in low b-value clusters. $M_{max}$ gives information on the largest rupture size during a burst, given that earthquake magnitude scales with the rupture area (Wells & Coppersmith, 1994). Similarly, the burst area can provide an upper bound for the largest rupture size of the events that occurred during the burst, assuming a burst represents the fault patch that activated following an episode of transient increase in shear stress. A scaling between the size of the activated area and the largest rupture size has been proposed for different types of seismic sequences, such as aftershock-mainshock sequences (Konstantinou et al., 2005) and those dominated by earthquake-to-earthquake interactions (Fischer & Hainzl, 2021).

Therefore, if $M_{max}$ and burst area remain small over time along a fault segment, this indicates that earthquake rupture size is inhibited, which is consistent with the associated high b-values, as this leads to a higher proportion of small events. On the other hand, if $M_{max}$ and burst area are large or increase over time, this indicates that earthquake rupture size is less constrained, which is consistent with the associated low b-values, as larger events can occur.

We point out that the activation of low b-value fault segments does not necessarily mean that large events will occur in the future. For example, cluster S1, one of the low b-value clusters, quickly deactivated in favor of high b-value faults that generated only low-magnitude seismicity (Figure 6). In addition, fault segments with similar extents can have significantly different largest magnitudes (e.g., N4 and S4, and S1 and S3; Figures 7c and 7d, and Figure S16 in Supporting Information S1). In other words, the total activated area is not always representative of the largest earthquake. Therefore, the rupture behavior of a fault segment, as inferred from the migration of seismicity, must be considered when predicting the largest rupture size and, consequently, assessing seismic hazard.

## 4.3. Physical Factors Affecting Event Rupture Size

To understand what physical processes may limit rupture size, we consider the main factors that are commonly invoked to explain variations in the b-value along different fault segments.

### 4.3.1. Stress State

The b-value of a fault system is typically interpreted to depend on its stress state. Many studies propose that the b-value is inversely proportional to stress, such as differential stress (Scholz, 1968, 2015) and shear stress (Mukuhira et al., 2024). A possible explanation for this relationship is that once a rupture starts propagating in a highly stressed region, it is less likely to be stopped by heterogeneities, which favors larger events and, in turn, lower b-values (e.g., Tormann et al., 2014).
Assuming that the eight fault segments are under the same strike-slip stress field, they are all well oriented for failure (Figures 4a and 4b). This suggests that the assumed stress field can provide a good first-order approximation of the stress acting on faults and can explain why seismicity tends to occur along two sets of orientations. However, this also means that the observed faults are under similar differential, normal, and shear stresses (Figure 4c). Therefore, the stress differences between clusters cannot explain the difference in the b-value. We note that, due to increased pore pressure from injection operations, clusters occurring during injection and near the injection points likely experience lower effective normal stress and, consequently, higher Coulomb stress than what is reported in Figure 4b. Nevertheless, clusters that occur after injection and farther from the injection points (e.g., southern clusters), where the pore pressure increase is expected to be low, still exhibit different b-values, in agreement with the conclusion that b-value variations do not result from differences in stress conditions along the fault segments. Therefore, other factors that can control the rupture size, and thus the b-value, must be considered.

### 4.3.2. Fault Properties

Changes in the b-value have also been attributed to variations in the structural and geomechanical properties of faults. For example, high-roughness fault systems can inhibit the formation of large events compared to smooth and mature faults, leading to high b-values (Cochran et al., 2023; Goebel et al., 2017). More generally, strength heterogeneity, and by extension, the spatial distribution of high-strength patches (also called strength barriers), has been proposed to influence rupture propagation and thus control the b-values (Steacy et al., 1996). In particular, faults with abundant or coherent high-strength patches will tend to have higher b-values than those where strength barriers are fewer in number, smaller, or weaker.

We point out that the presence of strength heterogeneity has also been proposed as a factor contributing to the discontinuous migration in time of the seismicity (Fischer et al., 2023), which is a key feature of the observed seismicity (Figure 8). We therefore propose that the presence of faults with different degrees of strength heterogeneity may be a key factor controlling fault rupture and, consequently, the b-value at the Geoven geothermal site.

The spatial proximity of low-b-value and high-b-value clusters indicates that significant changes in the fault geomechanical properties can occur over small distances. Moreover, the fact that all eight clusters are likely to be within the Hercynian granitic basement suggests that such variations can arise within the same lithology, highlighting the complex geology of the area.

## 5. Conclusions

Our analysis of the Vendenheim 2018–2022 sequence of induced seismicity revealed that seismicity primarily clustered along eight distinct fault segments. These faults ruptured intermittently through multiple bursts of seismicity, a behavior compatible with the triggering process of earthquake interactions.

The main observation of this study is that clusters occurring in close temporal and spatial proximity can have significantly different b-values. These variations can only be detected if the fault segments are properly identified through seismicity clustering analysis. This finding suggests that b-values are closely tied to the local physical conditions specific to each fault system. This conclusion is further supported by the fact that changes in the b-value are associated with different rupture behaviors along the faults: in high b-value faults, both burst area and $M_{max}$ remain small over time, whereas in low b-value faults, they can grow and become large.

The fault segments are well aligned with the regional stress field and, given their similar depths, are subject to comparable stress conditions. This indicates that while the regional stress field can explain why the activated faults have the observed orientation, it cannot account for variations in the b-values.

We propose that a possible factor that could have influenced rupture propagation, and consequently controlled the b-values along the fault segments, is the different degree of strength heterogeneity. High b-value clusters may correspond to fault segments with strength barriers (either determined by fault geometry or material properties) that inhibit rupture propagation, whereas low b-value clusters may represent fault segments where strength distribution is more homogeneous, which instead favors rupture propagation. The presence of strength barriers, including cohesion variability or complex hydro-thermal stresses, is consistent with the discontinuous migration of seismicity along the faults and the observation of smaller burst areas and $M_{max}$ in high b-value faults compared to low b-value faults.
Finally, this study emphasizes the importance of identifying and isolating the primary fault structures and their properties to better understand the spatiotemporal evolution of seismicity. More broadly, it highlights the value of high-resolution catalogs in advancing our understanding of fault rupture processes and their associated seismic hazards.

## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

In this study, we used waveform data from seismic networks with FDSN codes AM (Raspberry Shake, 2016), FR (Epos-France, 1962), GR (BGR, 1976), LE (Erdbebendienst Südwest & Rheinland-Pfalz, 2009), and RA (RSEIF, 1995). Waveform data from the Fonroche Géothermie stations, as well as from network XX, are available on Zenodo (Minetto, 2025b). All station codes are shown in Figure S2 of Supporting Information S1. The 3D GeORG geological model is accessible via the Cupidon platform (Capar, 2024). Version 1 of the DETX software was used for template matching (Lengliné, 2025). The software (version 1.0) for computing the stress acting on the faults is published on GitHub https://github.com/riccmin/fault-stress-calculator (Minetto, 2025a). The high-resolution seismicity catalog can be downloaded from the CDGP platform at https://cdgp.u-strasbg.fr/geonetwo-rk/srv/fre/catalog.search#/metadata/VDHD-Catalogue-Minetto (EOST & Fonroche Géothermie, 2020).

## References

Aki, K. (1965). Maximum likelihood estimate of b in the formula log N=a-bM and its confidence limits. *Bulletin of the Earthquake Research Institute (ERI) at the University of Tokyo*, **43**, 237–239. Retrieved from https://cir.nii.ac.jp/crid/1573387450038851840

Beeler, N. M., Simpson, R. W., Hickman, S. H., & Lockner, D. A. (2000). Pore fluid pressure, apparent friction, and Coulomb failure. *Journal of Geophysical Research*, **105**(B11), 25533–25542. https://doi.org/10.1029/2000JB900119

BGR. (1976). German Regional Seismic Network (GRSN) [Dataset]. Bundesanstalt für Geowissenschaften und Rohstoffe. https://doi.org/10.25928/mbx6-hr74

Byerlee, J. (1978). Friction of rocks. In J. D. Byerlee & M. Wyss (Eds.), *Rock friction and earthquake prediction* (pp. 615–626). Birkhäuser. https://doi.org/10.1007/978-3-0348-7182-2_4

Candela, T., Wassing, B., Ter Heege, J., & Buijze, L. (2018). How earthquakes are induced. *Science*, **360**(6389), 598–600. https://doi.org/10.1126/science.aat2776

Capar, L. (2024). Modèle géologique 3d de la partie française du fossé rhénan [Dataset]. BRGM. https://doi.org/10.18144/4a58a7e3-8b71-4c69-9447-cb1c791d23a3

Chen, H., Meng, X., Niu, F., Tang, Y., Yin, C., & Wu, F. (2018). Microseismic monitoring of stimulating shale gas reservoir in SW China: 2. Spatial clustering controlled by the preexisting faults and fractures. *Journal of Geophysical Research: Solid Earth*, **123**(2), 1659–1672. https://doi.org/10.1002/2017JB014491

Cochran, E. S., Page, M. T., van der Elst, N. J., Ross, Z. E., & Trugman, D. T. (2023). Fault roughness at seismogenic depths and links to earthquake behavior. *The Seismic Record*, **3**(1), 37–47. https://doi.org/10.1785/0320220043

Cochran, E. S., Ross, Z. E., Harrington, R. M., Dougherty, S. L., & Rubinstein, J. L. (2018). Induced earthquake families reveal distinctive evolutionary patterns near disposal Wells. *Journal of Geophysical Research: Solid Earth*, **123**(9), 8045–8055. https://doi.org/10.1029/2018JB016270

Collettini, C., & Tinti, E. (2025). The influence of lithology and fault source volume on the Magnitude-frequency-distribution of earthquakes. *Geophysical Research Letters*, **52**(2), e2024GL110354. https://doi.org/10.1029/2024GL110354

Cornet, F. H., Bérard, T., & Bourouis, S. (2007). How close to failure is a granite rock mass at a 5 km depth? *International Journal of Rock Mechanics and Mining Sciences*, **44**(1), 47–66. https://doi.org/10.1016/j.ijrmms.2006.04.008

De Barros, L., Cappa, F., Deschamps, A., & Dublanchet, P. (2020). Imbricated aseismic slip and fluid diffusion drive a seismic swarm in the Corinth Gulf, Greece. *Geophysical Research Letters*, **47**(9), e2020GL087142. https://doi.org/10.1029/2020GL087142

Dublanchet, P. (2020). Stress-dependent b value variations in a heterogeneous rate-and-state fault model. *Geophysical Research Letters*, **47**(13), e2020GL087434. https://doi.org/10.1029/2020GL087434

Efron, B., & Tibshirani, R. J. (1994). *An introduction to the bootstrap*. CRC Press. (Google-Books-ID: gLlplUxRntoC).

El-Isa, Z. H., & Eaton, D. W. (2014). Spatiotemporal variations in the b-value of earthquake magnitude-frequency distributions: Classification and causes. *Tectonophysics*, **615**, 1–11. https://doi.org/10.1016/j.tecto.2013.12.001

EOST, & Géothermie, F. (2020). Episode: Geoven geothermal project at vendenheim. Eost-cdgp [collection]. *EOST and Fonroche Géothermie*. https://doi.org/10.25577/kkz6-fc66

Epos-France. (1962). Epos-France broad-band network (rlbp) [Dataset]. *Epos-France Seismological Data Centre*. https://doi.org/10.15778/RESIF.FR

Ester, M., Kriegel, H.-P., Sander, J., & Xu, X. (1996). A density-based algorithm for discovering clusters in large spatial databases with noise. In *Proceedings of the second international conference on knowledge discovery and data mining* (pp. 226–231). AAAI Press.

Evans, K., Valley, B., Häring, M., Hopkirk, R. J., Baujard, C., Kohl, T., et al. (2009). *Studies and support for the EGS reservoirs at soultz-sous-forêts (technical report)*. Centre for Geothermal Research.

Fischer, T., & Hainzl, S. (2021). The growth of earthquake clusters. *Frontiers in Earth Science*, **9**, 638336. https://doi.org/10.3389/feart.2021.638336

Fischer, T., Hainzl, S., & Vlček, J. (2023). Fast migration episodes within earthquake swarms. *Geophysical Journal International*, **235**(1), 312–325. https://doi.org/10.1093/gji/ggad221

## Acknowledgments

We thank M. Grunberg, S. Lambotte, B. Derode, D. Zigone, F. Mattern, J. Vergne, A. Schlupp, H. Jund, F. Engels, S. Benlalam for very fruitful discussions. This research was funded in part or entirely by the Agence Nationale de la Recherche (ANR) under PreSENCE project ANR-21-CE05-0033. This work of the Interdisciplinary Thematic Institut Geosciences 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 SFRI-STRAT'US project (ANR-20-SFRI-0012).
Ge, S., & Saar, M. O. (2022). Review: Induced seismicity during geoenergy development—A hydromechanical perspective. *Journal of Geophysical Research: Solid Earth*, 127(3), e2021JB023141. https://doi.org/10.1029/2021JB023141

GFZGEMPA GmbH. (2008). The SeisComP seismological software package. GFZ Data Services. https://doi.org/10.5880/GFZ.2.4.2020.003

Gibbons, S. J., & Ringdal, F. (2006). The detection of low magnitude seismic events using array-based waveform correlation. *Geophysical Journal International*, 165(1), 149–166. https://doi.org/10.1111/j.1365-246X.2006.02865.x

Goebel, T. H., Kwiatek, G., Becker, T. W., Brodsky, E. E., & Dresen, G. (2017). What allows seismic events to grow big? Insights from b-value and fault roughness analysis in laboratory stick-slip experiments. *Geology*, 45(9), 815–818. https://doi.org/10.1130/G39147.1

Grigoli, F., Cesca, S., Priolo, E., Rinaldi, A. P., Clinton, J. F., Stabile, T. A., et al. (2017). Current challenges in monitoring, discrimination, and management of induced seismicity related to underground industrial activities: A European perspective. *Reviews of Geophysics*, 55(2), 310–340. https://doi.org/10.1002/2016R000542

Gulia, L., Rinaldi, A. P., Tormann, T., Vannucci, G., Enescu, B., & Wiemer, S. (2018). The effect of a mainshock on the size distribution of the aftershocks. *Geophysical Research Letters*, 45(24), 13277–13287. https://doi.org/10.1029/2018GL080619

Gutenberg, B., & Richter, C. F. (1944). Frequency of earthquakes in California. *Bulletin of the Seismological Society of America*, 34(4), 185–188. https://doi.org/10.1785/BSSA0340040185

Herrmann, M., Kraft, T., Tormann, T., Scarabello, L., & Wiemer, S. (2019). A consistent high-resolution catalog of induced seismicity in Basel based on matched filter detection and tailored post-processing. *Journal of Geophysical Research: Solid Earth*, 124(8), 8449–8477. https://doi.org/10.1029/2019JB017468

Herrmann, M., Piegari, E., & Marzocchi, W. (2022). Revealing the spatiotemporal complexity of the magnitude distribution and b-value during an earthquake sequence. *Nature Communications*, 13(1), 5087. https://doi.org/10.1038/s41467-022-32755-6

Igonin, N., Zecevic, M., & Eaton, D. W. (2018). Bilinear magnitude-frequency distributions and characteristic earthquakes during hydraulic fracturing. *Geophysical Research Letters*, 45(23), 12866–12874. https://doi.org/10.1029/2018GL079746

Jolliffe, I. T. (2002). *Principal component analysis*. Springer-Verlag. https://doi.org/10.1007/b98835

Keranen, K. M., & Weingarten, M. (2018). Induced seismicity. *Annual Review of Earth and Planetary Sciences*, 46(1), 149–174. https://doi.org/10.1146/annurev-earth-082517-010054

Kettlety, T., Verdon, J. P., Butcher, A., Hampson, M., & Craddock, L. (2020). High-resolution imaging of the ML 2.9 August 2019 earthquake in Lancashire, United Kingdom, induced by hydraulic fracturing during Preston new road PNR-2 operations. *Seismological Research Letters*, 92(1), 151–169. https://doi.org/10.1785/0220200187

Konstantinou, K. I., Papadopoulos, G. A., Fokaefs, A., & Orphanogiannaki, K. (2005). Empirical relationships between aftershock area dimensions and magnitude for earthquakes in the Mediterranean Sea region. *Tectonophysics*, 403(1), 95–115. https://doi.org/10.1016/j.tecto.2005.04.001

Lengliné, O. (2025). Detx [Software]. Zenodo. https://doi.org/10.5281/zenodo.15222295

Lengliné, O., Schmittbuhl, J., Drif, K., Lambotte, S., Grunberg, M., Kinscher, J., et al. (2023). The largest induced earthquakes during the Geoven deep geothermal project, Strasbourg, 2018–2022: From source parameters to intensity maps. *Geophysical Journal International*, 234(3), 2445–2457. https://doi.org/10.1093/gji/ggad255

Lomax, A., Virieux, J., Volant, P., & Berge-Thierry, C. (2000). Probabilistic earthquake location in 3D and layered models. In C. H. Thurber & N. Rabinowitz (Eds.), *Advances in seismic event location* (pp. 101–134). Springer Netherlands. https://doi.org/10.1007/978-94-015-9536-0-5

Minetto, R. (2025a). riccmin/fault-stress-calculator [Software]. Zenodo. https://doi.org/10.5281/zenodo.15222295

Minetto, R. (2025b). Vendenheim induced seismicity waveform dataset (2018–2022) [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.17111739

Mukuhira, Y., Fehler, M. C., Bjarkason, E. K., Ito, T., & Asanuma, H. (2024). On the b-value dependency of injection-induced seismicity on geomechanical parameters. *International Journal of Rock Mechanics and Mining Sciences*, 174, 105631. https://doi.org/10.1016/j.ijrmms.2023.105631

Muntendam-Bos, A. G., & Grobbe, N. (2022). Data-driven spatiotemporal assessment of the event-size distribution of the Groningen extraction-induced seismicity catalogue. *Scientific Reports*, 12(1), 10119. https://doi.org/10.1038/s41598-022-14451-z

Norris, J. Q., Turcotte, D. L., Moores, E. M., Brodsky, E. E., & Rundle, J. B. (2016). Fracking in tight shales: What is It, what does it accomplish, and what are its consequences? *Annual Review of Earth and Planetary Sciences*, 44(1), 321–351. https://doi.org/10.1146/annurev-earth-060115-012537

Ogata, Y., & Katsura, K. (1993). Analysis of temporal and spatial heterogeneity of magnitude frequency distribution inferred from earthquake catalogues. *Geophysical Journal International*, 113(3), 727–738. https://doi.org/10.1111/j.1365-246X.1993.tb04663.x

Raspberry Shake, S. (2016). Raspberry shake [Dataset]. *International Federation of Digital Seismograph Networks*. https://doi.org/10.7914/SN/AM

Reynolds, D. (2009). Gaussian mixture models. In S. Z. Li & A. Jain (Eds.), *Encyclopedia of biometrics* (pp. 659–663). Springer US. https://doi.org/10.1007/978-0-387-73003-5_196

Rothé, J.-P., & Peterschmitt, E. (1950). Etude séismique des explosions d'Haslach. *Annales de l'Institut de physique du Globe de Strasbourg*, 5(1), 13–38. https://doi.org/10.3406/aipgs.1950.1881

RSEIF. (1995). Resif-rap French accelerometric network [Dataset]. RESIF—Réseau Sismologique et géodésique Français. https://doi.org/10.15778/RESIF.RA

Sáez, A., & Lecampion, B. (2023). Post-injection aseismic slip as a mechanism for the delayed triggering of seismicity. *Proceedings of the royal society A: Mathematical, physical and engineering sciences* (Vol. 479(2273), p. 20220810). https://doi.org/10.1098/rspa.2022.0810

Schmittbuhl, J., Lambotte, S., Lengliné, O., Grunberg, M., Jund, H., Vergne, J., et al. (2021). Induced and triggered seismicity below the city of Strasbourg, France from November 2019 to January 2021. *Comptes Rendus. Geoscience*, 353(S1), 561–584. https://doi.org/10.5802/crgeos.71

Scholz, C. H. (1968). The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes. *Bulletin of the Seismological Society of America*, 58(1), 399–415. https://doi.org/10.1785/BSSA0580010399

Scholz, C. H. (2015). On the stress dependence of the earthquake b value. *Geophysical Research Letters*, 42(5), 1399–1402. https://doi.org/10.1002/2014GL062863

Shapiro, S. A., Huenges, E., & Borm, G. (1997). Estimating the crust permeability from fluid-injection-induced seismic emission at the KTB site. *Geophysical Journal International*, 131(2), F15–F18. https://doi.org/10.1111/j.1365-246X.1997.tb01215.x

Shelly, D. R., Ellsworth, W. L., & Hill, D. P. (2016). Fluid-faulting evolution in high definition: Connecting fault structure and frequency-magnitude variations during the 2014 long valley Caldera, California, earthquake swarm. *Journal of Geophysical Research: Solid Earth*, 121(3), 1776–1795. https://doi.org/10.1002/2015JB012719

Spada, M., Tormann, T., Wiemer, S., & Enescu, B. (2013). Generic dependence of the frequency-size distribution of earthquakes on depth and its relation to the strength profile of the crust. *Geophysical Research Letters*, 40(4), 709–714. https://doi.org/10.1029/2012GL054198
Steacy, S. J., McCloskey, J., Bean, C. J., & Ren, J. (1996). Heterogeneity in a self-organized critical earthquake model. *Geophysical Research Letters*, 23(4), 383-386. https://doi.org/10.1029/96GL00257

Stormo, A., Lengliné, O., & Schmittbuhl, J. (2015). Mechanical origin of b-value changes during stimulation of deep geothermal reservoirs. *Geothermal Energy*, 3, 1-18. https://doi.org/10.1186/s40517-014-0022-0

Südwest, E., & Rheinland-Pfalz. (2009). Erdbebendienst südwest [Dataset]. *International Federation of Digital Seismograph Networks*. https://doi.org/10.7914/SN/LE

Terrier, M., de Santis, F., Soliva, R., Valley, B., Bruel, D., Geraud, Y., & Schmittbuhl, J. (2022). Rapport Phase 1 du comité d'experts créé en appui à l'administration sur la boucle géothermique GEOVEN (Tech. Rep.).

Tormann, T., Wiemer, S., & Mignan, A. (2014). Systematic survey of high-resolution b value imaging along Californian faults: Inference on asperities. *Journal of Geophysical Research: Solid Earth*, 119(3), 2029-2054. https://doi.org/10.1002/2013JB010867

Trugman, D. T., & Shearer, P. M. (2017). GrowClust: A hierarchical clustering algorithm for relative earthquake relocation, with application to the Spanish springs and Sheldon, Nevada, earthquake sequences. *Seismological Research Letters*, 88(2A), 379-391. https://doi.org/10.1785/0220160188

van der Elst, N. J. (2021). B-Positive: A robust estimator of aftershock magnitude distribution in transiently incomplete catalogs. *Journal of Geophysical Research: Solid Earth*, 126(2), e2020JB021027. https://doi.org/10.1029/2020JB021027

Wells, D. L., & Coppersmith, K. J. (1994). New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. *Bulletin of the Seismological Society of America*, 84(4), 974-1002. https://doi.org/10.1785/BSSA0840040974

Woessner, J., & Wiemer, S. (2005). Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty. *Bulletin of the Seismological Society of America*, 95(2), 684-698. https://doi.org/10.1785/0120040007

Wyss, M., et al. (1973). Towards a physical understanding of the earthquake frequency distribution. *Geophys. JR Astron. Soc*, 31(4), 341-359. https://doi.org/10.1111/j.1365-246x.1973.tb06506.x

Wyss, M., Hasegawa, A., Wiemer, S., & Umino, N. (1999). Quantitative mapping of precursory seismic quiescence before the 1989 M 7.1 off-Sanriku earthquake, Japan. *Annals of Geophysics*, 42(5), 20. https://doi.org/10.4401/ag-3765

Zoback, M. D. (2007). *Reservoir geomechanics*. (pp. 154-157). Cambridge University Press. https://doi.org/10.1017/CBO9780511586477