# The largest induced earthquakes during the GEOVEN deep geothermal project, Strasbourg, 2018-2022: from source parameters to intensity maps

Olivier Lengliné, Jean Schmittbuhl, Karim Drif, Sophie Lambotte, Marc Grunberg, Jannes L. Kinscher, Christophe Sira, Antoine Schlupp, Marc Schaming, Hélène Jund, et al.

► To cite this version:

Olivier Lengliné, Jean Schmittbuhl, Karim Drif, Sophie Lambotte, Marc Grunberg, et al. The largest induced earthquakes during the GEOVEN deep geothermal project, Strasbourg, 2018-2022: from source parameters to intensity maps. Geophysical Journal International, 2023, 234, pp.2446-2458. 10.1093/gji/ggad255. insu-04198297

HAL Id: insu-04198297

https://insu.hal.science/insu-04198297

Submitted on 7 Sep 2023

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

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

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

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

# The largest induced earthquakes during the GEOVEN deep geothermal project, Strasbourg, 2018–2022: from source parameters to intensity maps

O. Lengliné<sup>1</sup>, J. Schmittbuhl,<sup>1</sup> K. Drif,<sup>1,2</sup> S. Lambotte,<sup>1</sup> M. Grunberg,<sup>3</sup> J. Kinscher,<sup>2</sup> C. Sira,<sup>3</sup> A. Schlupp,<sup>1</sup> M. Schaming,<sup>1</sup> H. Jund<sup>3</sup> and F. Masson<sup>1</sup>

Accepted 2023 June 22. Received 2023 June 19; in original form 2022 December 5

## SUMMARY

Between 2019 November and 2021 July, four induced earthquakes of local magnitude equal to or greater than three were felt by the population of Strasbourg, France. These events were related to activity at the deep geothermal site GEOVEN located in Vendenheim in the northern suburb area of the city of Strasbourg. The first earthquake, with a local magnitude ($M_{Iv}$) of 3.0, occurred on 2019 November 12, at the same depth as the bottom of the wells (approximately 4 km) but 5 km to the south. The second ($M_{Iv}$ 3.6) occurred a year later, on 2020 December 4, below the wells, and led to the termination of the project by the authorities. The third ($M_{Iv}$ 3.3) was initiated three weeks after shut-in on 2021 January 22, while the largest earthquake to date ($M_{Iv}$ 3.9) occurred on 2021 June 26, more than 6 months after shut-in. We constrained these four events' absolute locations using a 3-D velocity model of the area and here present regional intensity maps. We estimated moment magnitude and focal mechanism trough waveform inversion and inferred the fault plane activated during the largest event from an analysis of rupture directivity effects in the recorded waveforms. Our analysis highlights the existence of a critically stressed fault that hosted three of these widely felt events. We show how the derived source properties of these four earthquakes are directly linked to ground shaking observations at the surface. Notably, we demonstrate how earthquake moment, location, direction of rupture and stress drop impact the regional intensity distribution. Our results suggest that the traffic light system could benefit from including ground shaking scenarios based on realistic subsurface properties and potential earthquake source models.

**Key words:** Waveform inversion; Earthquake source observations; Induced seismicity.

## 1 INTRODUCTION

A series of more than 500 earthquakes occurred below the northern suburb of the city of Strasbourg between the end of 2019 and 2021 July (Schmittbuhl et al. 2021). They were recorded and localized by the French service for the observation of national seismic activity, hosted by EOST (BCSF-Rénass). Four of these events were particularly felt by the local population. The first event, with a local magnitude of $M_{Iv}$ 3.0, occurred on the 2019 November 12 and was located below the Robertsau area. The second event occurred on the 2020 December 4 with a local magnitude of $M_{Iv}$ 3.6. It was located 5.5 km further north under the city of La Wantzenau, close to the termination of one of the deep wells of the deep geothermal energy project GEOVEN. The third event was close to the $M_{Iv}$ 3.6

event, occurring on the 2021 January 21, with a local magnitude of $M_{Iv}$ 3.3. The last, the largest of the sequence thus far, took place close to the two previous events and had a magnitude $M_{Iv}$ 3.9.

These four events were unexpected. No natural earthquakes had been recorded over the previous 40 yr within a 5 km radius around the GEOVEN project before their onset (Schmittbuhl et al. 2021). Given that the area's natural seismicity is very low and that all these earthquakes occurred after the onset of the GEOVEN project and in close vicinity to the boreholes in which fluid was injected, these events are illustrative examples of induced felt earthquakes in an urban context that strongly impacted the local population. The events were much greater than the alert level of the traffic light system (TLS) implemented during the GEOVEN project. Industrial operators and regulators are typically looking for tools and procedures to

<footnote n="1"><sup>1</sup>Université de Strasbourg, CNRS, ITES, UMR7063, 5 rue René Descartes, 67084 Strasbourg Cedex, France. E-mail: lengline@unistra.fr</footnote>

<footnote n="2"><sup>2</sup>INERIS, Campus ARTEM, 92 rue de Sergent Blandan, 54042 Nancy Cedex, France</footnote>

<footnote n="3"><sup>3</sup>Université de Strasbourg, CNRS, EOST/UAR 830, 5 rue René Descartes, 67084 Strasbourg Cedex, France</footnote>
control seismicity and assess the related hazards, using TLS systems with typically three levels (green—yellow—red) corresponding to increasing thresholds of potential vibrations (Bommer & Crowley 2006; Grigoli et al. 2017; Mignan et al. 2017). TLS can be physics-based (e.g. Kwiatek et al. 2019) but are generally based on simple parameters such as magnitude, maximum peak ground acceleration (PGA) or maximum peak ground velocity (PGV) according to the emphasis given to the hazard view, that is, the characterization of the seismic source in the reservoir (magnitude of the event) or the response of the ground below buildings (PGV or PGA). Here, the red level of the TLS was set using a magnitude criterion of $M_{lv} = 2$.

Despite the seismic hazard posed by induced earthquakes (Foulinger et al. 2018), these events are useful for assessing the development of the reservoir, as they can delineate the activated structures (e.g. Lengliné et al. 2017) in the reservoir, provide key indicators of the underground state of stress (Cornet & Jianmin 1995; Schoenball et al. 2014), track the evolution of the pore pressure (e.g. Shapiro et al. 2002), or help to infer geomechanical parameters of the reservoir (e.g. McClure 2012).

Whether the objective is gaining insight into reservoir behaviour or implementing a TLS, both approaches rely on the same measurements: the recorded waveforms of monitoring seismic stations. However, the design of a monitoring seismic network might be different for each approach: as sensitive as possible for describing small events and for reservoir management (dense network, borehole arrays, low noise site, highly sensitive sensors, etc.), or as close as possible to infrastructure for site effect assessment (Zang et al. 2014; Bohnhoff et al. 2018). In practice, monitoring networks are typically a compromise for assessing both objectives, and TLSs are evolving to be adaptive, including updated on-the-fly geomechanical and seismological parameters and, ideally, attempting to include the transfer function of the dynamic strain from the seismic source to the ground response (Grigoli et al. 2017; Mignan et al. 2017). However, advanced TLS systems (able to predict seismicity and account for variable scenarios) are limited in their ability to include multiple records, advanced analysis and detailed indicators in quasi real time. Moreover, TLSs rely on the principle that seismicity is not only important at the origin of the risk when events are large, but is also, the key indicator when events are small, since seismicity can provide in-depth knowledge of the dynamics of the reservoirs and can potentially forecast the nucleation phase of the largest events. This principle relies on numerous assumptions, including: high sensitivity and reliability in monitoring systems, appropriate description of the velocity model around the reservoir and a seismic signature of the main deformation of the reservoir. These assumptions are increasingly questioned with respect to important aseismic deformation (Cornet et al. 1997; Schmittbuhl et al. 2014; Wei et al. 2015; Guglielmi et al. 2015; Cornet 2016; Hopp et al. 2019). This work aims to improve the link between ground shaking at the surface and seismic sources at depth. A better knowledge of this link will have strong implications for improving reservoir monitoring and optimizing TLS strategies.

Our first objective is to decipher the main features of the four largest induced earthquakes that occurred in the Strasbourg sequence. We assess their characteristics in terms of surface observations, such as shaking intensity maps and their seismic source properties at depth. We then aim to link these features by developing a 3-D velocity model around the reservoir with deterministic waveform modelling. Thus, our analysis demonstrates how earthquake sources impacts on the population could be calculated and integrated into TLS.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Representation of S-wave velocity profiles at three locations in the Rhine Graben. All depths are relative to sea level. The green dashed line shows the 1-D reference velocity (similar for all locations). The blue line corresponds to velocity profiles extracted from our 3-D model. The red line gives the reference velocity model derived at Soultz-sous-Forêts (Charléty et al. 2006) and Rittershoffen (Maurer et al. 2020). The yellow and grey areas represent the sedimentary cover and the granitic basement, respectively.</figcaption></figure>

## 2 EARTHQUAKE LOCATIONS

### 2.1 3-D velocity model of area

A suitable velocity model for the reservoir area and its surroundings is required to obtain a precise estimate of the event locations and, subsequently, of the seismic source properties of the events. The BCSF-Rénass's default velocity model for the region is a 1-D velocity model based on Rothé & Peterschmitt (1950), which does not accurately represent the shallow sedimentary layers of the graben. Sedimentary layers can be quite thick in the upper Rhine graben, particularly in the city of Strasbourg and its surroundings. A set of regional normal faults with variable offsets also contribute to significant lateral variability in the basin. As a result, there can be significant velocity contrasts between these layers and the 1-D reference velocity model which could result in complex patterns on the intensity maps. Here, we created a 3-D velocity model by considering the varying thickness of the Rhine graben's sediment cover. Based on the data set released by Freymark et al. (2020), which maps all geological units in the upper Rhine Graben with their appropriate depth and thickness, the 3-D model has a variable basement depth. In particular, we identified the transition depth and mapped (on a horizontal grid with a resolution of 1 km) the bottom of the Buntsandstein unit of this data set, which represents the interface between the sedimentary cover and the granitic basement. We considered that the regional 1-D reference velocity model (Fig. 1) could accurately describe all seismic velocities below this transition depth. We decided to set the *P*- and *S*-wave velocities in the sedimentary cover, above the transition depth, as varying linearly in steps of 500 m from the surface to the transition depth, with velocities near the surface fixed at 2.0 and 0.9 km s⁻¹ and those at the top of the granitic basement at 5.2 and 3.1 km s⁻¹ for *P*
and S waves, respectively. These numbers are extrapolated from the velocity models at the two geothermal sites in Rittershoffen and Soultz-sous-Forêts. A comparison of our model and the reference 1-D velocity model shows some notable differences (Fig. 1). We also note that the two velocity models of Rittershoffen (Maurer et al. 2020) and Soultz-sous-Forêts (Charléty et al. 2006) are well captured by our straightforward linear gradient 3-D velocity model. In fact, 1-D profiles at the locations of these two reservoirs extracted from our 3-D model agree well with the measured profile at each site. Due to the orientation of the Rhine graben, the 3-D model exhibits an N-S direction with slow velocity, as shown in Fig. S1 (Supporting Information). We also extract a 1-D profile at the location of the GEOVEN wells (Vendenheim profile in Fig. 1), and we note that the transition between the granitic basement and the sedimentary cover is, on average, deeper here than at the Soultz-sous-Forêts and Rittershoffen sites, at 3.5 km depth, indicating the presence of a thickening of the sedimentary cover at the Strasbourg site.

## 2.2 Earthquake locations in the 3-D model

When calculated using the 1-D reference model, the locations of the four largest events obtained by the BCSF-Rénass are related to some significant traveltime residuals at several stations, particularly for S-wave arrivals (Fig. S2, Supporting Information). These residuals point to a significant difference in seismic wave velocity between the various stations, which may have an impact on the earthquakes’ inverted locations. Due to the heterogeneous nature of the subsurface, we use the 3-D velocity model that was developed in the previous section to both better constrain the earthquake locations and minimize these effects. Using NonLinLoc software (Lomax et al. 2001), we invert the locations of the events based on the manual *P* and *S* wave arrival time picks made by BCSF-Rénass operators. In Fig. 2, we show the location and related location uncertainty for each of these events. When compared to the results obtained using the reference 1-D velocity model, the traveltime residuals from our location in our 3-D model are significantly reduced (Fig. S2, Supporting Information). In fact, we reduce the mean *S*-wave residual from 0.57 to 0.17 s and the mean *P*-wave residual from 0.23 to 0.10 s. These four events all have depths of approximately 4300 m, which is roughly equal to the depth of the wells’ bottoms but is 1 km shallower than the BCSF-Rénass locations (Fig. S3, Supporting Information). Our findings also show that the three most recent events occurred close to the GEOVEN wells, specifically between the two wells. Our locations move these earthquakes 1 km south of where the BCSF-Rénass places them. These three events are aligned in an NW-SE direction. We note that this direction is resolved in Schmittbuhl et al. (2021) only from the relative relocation of earthquakes.

# 3 INTENSITY MAPS

## 3.1 Macroseismic and instrumental data

For all four earthquakes, spontaneous testimonies from the public arrived in the first minutes on the www.franceseisme.fr website. An extensive call for testimonies was made through social networks, the media and the BCSF-Rénass seismological information distribution list. In addition to questions regarding the impacts on common indicators (people, objects, furnitureandildings), the individual forms filled out by the witnesses, were accompanied by representative and

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Location of the four main events based on picked *P*- and *S*-wave arrival times and using the built 3-D velocity model (stars). The ellipses around each star represent the location uncertainties. The colour refers to each earthquake: red is the 2019 November event, blue is the 2020 December event, green is the 2021 January event and yellow is the 2021 June event. The plain circles refer to the location of the same events by the BCSF-Rénass. The brown lines show the simplified paths of the geothermal wells. The purple lines indicate the known fault location at the top of the basement in the GEORG model. The blue triangles show the locations of the seismic stations used in this study.</figcaption></figure>

selectable thumbnails of seismic effects that provided an individual value of ground shaking severity (Sira et al. 2021b). Averaged over the municipality, they gave a rapid and preliminary measure of intensity that is fairly reliable for intensities below V (strong shaking) when the number of returns is sufficiently large (over 10 forms collected, Sira et al. 2021a).

Except for the earthquake of 2019 November 12, for which only these individual forms were used, BCSF-Rénass also took into account, in its estimation of final municipal intensities, the results of collective forms returned by town halls. For these forms, the BCSF-Rénass has generated a specific building indicator for reporting damage to buildings by inhabitants. Only the municipality of La Wantzenau (1 km to the epicentres of the 2020 and 2021 earthquakes) returned a significant number of declared damages to the BCSF-Rénass. The intensity estimates were made according to the criteria of the EMS-98 macroseismic scale (Grünthal 1998) essential for estimating the severity of shaking.

For German territory, the Landesamt für Geologie, Rohstoffe und Bergbau (LGRB) and the Landeserdbebendienst of Baden-Württemberg, and the Landesamt für Geologie und Bergbau (LGB) of Rhineland-Palatine also collected the testimonies of the inhabitants through their websites. Based on these testimonies, we estimated an EMS-98 intensity for localities with more than 10 collected individual forms. These intensities have been added to the cross-border macroseismic maps produced for the four main Strasbourg earthquakes (Fig. 3 and Table 1, and Fig. S4, Supporting Information).

The three macroseismic studies (for the 2020 December 04, 2021 January 22 and 2021 June 26 earthquakes) indicate maximum intensities of IV–V to V for the French and German territory near the epicentres. The statistical robustness of the EMS98 scale suggests
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Intensity maps for the four studied events. For the 2019 event, the intensity map is from citizen testimonies received by the BCSF-Rénass and LGRB, averaged by localities. For the three others, intensities are estimated using communal and individual forms for the French portion and only individual forms from localities in Germany. The red stars indicate the locations of the epicentres.</figcaption></figure>

<table><caption>**Table 1.** Number of reports collected in the cross-border macroseismic survey for each of the four analysed earthquakes. The reports are listed according to their geographical origin (B.W. is Baden-Württemberg and R.P. Rhineland-Palatine).</caption><thead><tr><td>Date</td><td>Magnitude</td><td>France</td><td>B.W.</td><td>R.P.</td></tr></thead><tbody><tr><td>2019/11/12</td><td>M<sub>IV</sub> 3.0</td><td>555</td><td>133</td><td>0</td></tr><tr><td>2020/12/04</td><td>M<sub>IV</sub> 3.6</td><td>1540</td><td>977</td><td>0</td></tr><tr><td>2021/01/14</td><td>M<sub>IV</sub> 3.3</td><td>1014</td><td>1459</td><td>0</td></tr><tr><td>2021/06/26</td><td>M<sub>IV</sub> 3.9</td><td>1861</td><td>1450</td><td>3</td></tr></tbody></table>

with good certainty that the V intensity was not exceeded for any of the earthquakes. Most of the damage was of grade 1 on buildings with a vulnerability of A to C. Grade 2 occurred very rarely (<1 per cent of the buildings in each vulnerability class B and C). Combining this observation with other indicators is in agreement with an upper bound estimation of the intensity to V. The estimation must consider the coherence of all the indicators. Particularly in the case of an induced earthquake where the residents are well aware of a potential reimbursement, an excessive reliance on declared damages as a diagnosis results in an overestimation of intensity values. The damage data collected by the BCSF-Rénass after the earthquake
of 2021 June 26 are rather scarce and probably do not represent all effects produced on buildings (Sira & Grunberg 2021).

In addition to the intensity estimated from macroseismic studies and individual testimonies, we computed the intensity derived from PGV at instruments in the vicinity of the epicentre up to a distance of 60 km. We estimated PGV at each station using the maximum velocity amplitude on horizontal components when available. We used the relationship provided by Caprio et al. (2015) to convert these PGV values into intensities.

## 3.2 Intensity distribution

Generally, for the same magnitude, shallow earthquakes generate felt vibrations at shorter distances from the epicentre than deep earthquakes. For the analysed Strasbourg events, a notable vibration was produced up to distances of 120 km in both France and Germany, particularly for the earthquake of 2021 June 26 (Fig. S4, Supporting Information). The location of the witnesses in single-family homes or ground floors at these farthest distances suggests potential geological or topographical site effects, possibly related to the Rhine Graben's sedimentary terrain.

We compared the decay of the intensities with distance from the epicentre between each earthquake to reveal their respective impacts. We observe that there is good agreement between the intensity estimated from the macroseismic data and that derived from the PGV (Fig. 4). The largest magnitude earthquake of the sequence, the 2021 June, $M_{lv}$ 3.9 event, produced the greatest intensities at the farthest distances, as expected. The intensities linked to the 2020 December and 2021 January earthquakes are nearly identical and only slightly lower than those of the magnitude $M_{lv}$ 3.9 event. The lowest intensities are observed for the 2019 November earthquake. For this event, although its estimated magnitude is not very different from that of the 2021 January earthquake, the differences in intensities are quite pronounced.

# 4 EARTHQUAKE SOURCE PROPERTIES

## 4.1 Moment tensor inversion

Establishing the source properties of the four earthquakes is not a trivial task since they are small events (i.e. small magnitude) and occurred in an urban environment, which creates significant noise on the seismograms. Nevertheless, there exists a dense coverage of seismic instruments, of various qualities close to the epicentres (Fig. 2 and Fig. S5, Supporting Information). We estimated the focal mechanisms for each of the four events by inverting the components of their moment tensor based on the displacement waveforms recorded in the vicinity of the epicentre. By successively setting each of the components of the moment tensor to unity, we calculated the elastodynamic wave propagation in a 3-D medium for each source location using the finite-difference method of Maeda et al. (2017) (Text S1 and Table S1, Supporting Information). The velocity model established in the preceding section is used in the calculation. We used a grid search approach to invert the components of the moment tensors of these four main earthquakes by varying the strike, dip and rake, each in 1° steps, restricting our analysis to the search for double-couple mechanisms. We obtained the three-component synthetic displacement seismograms at each chosen site by summing the computed signals for each component of the moment tensor in accordance with Aki & Richards (2002). The displacement seismograms are filtered between 0.14 and 0.8 Hz and

with a duration of 10 s starting at the origin time of the earthquake. For the 2021 June 26 earthquake, we extend the duration to 15 s because of the longer duration of the signal. The frequency range and the signal duration are chosen based on the parameter settings from Wang & Zhan (2020) for earthquakes of similar magnitudes. The noise level at low frequencies determines the lowest frequency range, and the point source approximation and unmodelled small-scale subsurface heterogeneities determine the highest bound. The time window includes both the P and S waves in the same window (Fig. 5) due to the proximity of the earthquakes to the recording sites. By applying a time-shift to all traces that best maximizes the average correlation coefficient between synthetic and observed signals, we allowed for an uncertainty of the origin time of each event (Table S2, Supporting Information). We estimated the moment of each event from the difference in amplitude between the synthetics and the observed waveforms for each trial in the grid search. Finally, to assess the quality of the source parameters, we computed an L2 norm between the moment-scaled synthetics and the recorded signals after correcting for time-shift. We kept the earthquake locations determined using the 3-D velocity model. Given the small size of the events (as hypothesized from their magnitude), the difference in location between the centroid and hypocentre should be small, and this approximation is reasonable.

We report the best models for each event as defined as those for which the misfit is no more than 5 per cent of the best-fitting model (Fig. 6). The range of moment magnitudes associated with each of these best mechanisms is also determined (Table 2). We also report an independent validation of the obtained result for each event by plotting the P-wave polarities determined manually by BCSF-Rénass analysts. We found good agreement between these polarities and the recovered mechanisms by the waveform fitting scheme.

## 4.2 Spectral properties of the main events

We present an analysis of the spectra of the four studied events to determine the source size of the events and validate our moment estimates. For each event, at each station, a window of five s begins one s before the S-wave pick is extracted. For stations with no S-wave pick available, a theoretical S-wave arrival is computed in the 3-D velocity model based on the source and station location. The displacement spectrum, $U(f)$, is obtained after deconvolving from the instrument response and filtered in the 0.1–35 Hz frequency range (Fig. 7). For each station, when three component records are available in the spectrum, $U(f)$ is computed as $U(f) = \sqrt{U_z^2(f) + U_n^2(f) + U_c^2(f)}$, where $U_x(f)$ is the S-wave spectrum for component x. The spectrum at each station is then fitted in the [0.3,30] Hz range using Brune's (1970) spectral model and considering the attenuation, corresponding to

$$U^{th}(f) = \frac{\Omega_0}{1 + \left(\frac{f}{f_c}\right)^2} \exp(-\pi f t^*), \quad (1)$$

where $f_c$ is the corner frequency, $t^* = \frac{t}{Q_0}$ with $t$ being the S-wave traveltime between the source and the station, and $Q_0$ is the S-wave quality factor. The seismic moment, $M_0$, is deduced from $\Omega_0$ (Text S2, Supporting Information). The inversion uses a truncated Newton algorithm with bounds and minimizes the least-squares norm between the calculated and observed spectra (Satriano 2021). A weighted average estimate of each of the three parameters obtained at all stations is returned. We observe that our moment estimates are in close agreement with those inferred from the waveform fitting
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Macroseismic intensities (circles) and computed intensities derived from PGV (triangles) for each of the four studied earthquakes. For each event, we compute the average distance for each intensity level reported by macroseismic data (dashed line) and its median absolute deviation (coloured area). To better visualize the differences in intensities among events, we report in the upper left panel an insert where we plot this average tendency for each earthquake.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of the recorded displacement waveforms for selected stations of the earthquake of 2021 June (grey line). The portion of each waveform used for the inversion is shown in black. The signals are filtered between 0.14 and 0.8 Hz. The red curves show the synthetic waveforms obtained for the best-fitting moment tensor. For each waveform segment, the correlation coefficient between the synthetic and the observed waveform is reported along with the station name and the component.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Focal mechanisms for the four largest events. Colours are similar to Fig. 2. For each event, the best solutions are displayed, and P-wave polarities are represented. The purple nodal planes are those of the best-fitting solution. Purple circles surrounding polarities, and purple crosses indicate stations used for the inversion, with or without a determined polarity, respectively. The date of each event is indicated.</figcaption></figure>

<table><caption>**Table 2.** Table summarizing the results of the source parameters for the four analysed earthquakes. For each earthquake we indicate the moment magnitude, as inferred from the inverted moment tensor, $M_w$ or from spectral fitting, $M_w^*$. The magnitude has been calculated from the moment using the relation of Hanks & Kanamori (1979). The two nodal planes that are associated with the mechanism with the lowest misfit are called NP1 and NP2 and are described by their strike/dip/rake. We report as well the corner frequency, $f_c$ as deduced from spectral fitting.</caption><thead><tr><th>Date</th><th>M<sub>w</sub><sup>min</sup> - M<sub>w</sub><sup>max</sup></th><th>NP1</th><th>NP2</th><th>M<sub>w</sub><sup>*</sup></th><th>f<sub>c</sub> (Hz)</th></tr></thead><tbody><tr><td>2019/12/11</td><td>2.7-2.9</td><td>202/56/-16</td><td>301/77/-144</td><td>2.8 ± 0.3</td><td>5 ± 2</td></tr><tr><td>2020/12/04</td><td>3.1-3.2</td><td>36/85/8</td><td>305/82/175</td><td>3.1 ± 0.2</td><td>7 ± 2</td></tr><tr><td>2021/01/22</td><td>3.0-3.1</td><td>212/80/-20</td><td>306/70/-169</td><td>3.0 ± 0.3</td><td>5 ± 2</td></tr><tr><td>2021/06/26</td><td>3.4-3.5</td><td>124/88/-169</td><td>33/79/-2</td><td>3.4 ± 0.3</td><td>5 ± 2</td></tr></tbody></table>

and that the corner frequencies are quite similar for the four events (Table 2).

## 4.3 Directivity analysis of the 2021 June earthquake rupture

Schmittbuhl et al. (2021) identified a regional fault plane, included in the public database GEORG (http://www.geopotenziale.org) connecting the 2019 southern event and the northern cluster close to the wells, the so-called Robertsau fault. The direction of this fault varies between N025°E and N010°E from south to north (Fig. 2). This direction might match one of the nodal planes of the earthquake that occurred in 2019 November (Fig. 6). However, the seismicity distribution near the injection wells, as determined by double-difference relocation, reveals a preferred direction that mainly favours the auxiliary plane with N304°E (Schmittbuhl et al. 2021).

A geomechanical stability analysis of the faults can benefit from knowing which of these two nodal planes is the real fault plane.

Looking at potential azimuthal variation in the seismic radiation, or directivity effects of the rupture (Haskell 1964), is one way to demonstrate which of these two planes is the most appropriate. We examine the variation in the apparent rupture duration at the various recording sites to determine the rupture direction during the largest event of the sequence, the 2021 June event. We anticipate that, in the case of a unilateral rupture, the apparent rupture duration will vary depending on the angle formed by the rupture direction and ray direction. The stations in the direction of the rupture will have a lower apparent rupture duration, while in the opposite direction to the direction of rupture, we expect a higher apparent rupture duration.

To retrieve apparent rupture durations, we first calculate the apparent source time functions using the projected Landweber deconvolution method, which is effective in producing stable results (Bertero et al. 1997; Lanza et al. 1999; Vallée 2004). This approach relies on the deconvolution of the main shock signal by an approximation of the Green function of the medium taken as the signal of a colocated earthquake of smaller size. We found a suitable event
<figure><img src="image_9.png" /><figcaption>**Figure 7.** S-wave source spectra for the event of 2021 June recorded at Station R4B4A (the closest to hypocentre) of the three components (coloured plain lines) and noise spectra (dashed lines). The best-fitting Brune's model is displayed as a black line, while we also show the best fit without considering attenuation (grey solid line) or without corner frequency (grey dashed line). The vertical line shows the estimated corner frequency with its associated uncertainty (grey area). The horizontal grey area shows the estimated range of the moment.</figcaption></figure>

that can be used as an empirical Green's function (EGF) that occurred on 2021 April 10, with magnitude $M_{Iv}$ 2.7. The projected Landweber method is an iterative process that computes an estimate of the apparent source time function (astf) signal, $f(t)$. The projection enforces a positivity constraint as well as a constraint of a finite rupture duration, $T$, on the recovered astf signal (Text S3 and Fig. S6, Supporting Information).

To interpret the results of the deconvolution in terms of rupture direction, we need to compute the apparent rupture duration at each recording site. We estimate these apparent rupture durations by defining the onset and the end of the rupture as the points where the moment, as computed from the time integration of the astf, becomes greater than 5 and 95 per cent of the final moment, respectively. The apparent rupture duration, $\tau_i$, at station $i$ is finally obtained as the duration between these two times.

According to the Haskell (1964) rupture model and following Park & Ishii (2015) and Abercrombie et al. (2017), the theoretical apparent rupture duration, $\tau_i^{th}$, recorded at station $i$ will vary with

$$ \tau_i^{th} = \tau_0 + \frac{L}{V_i} [\sin(\gamma_r) \sin(\gamma_i) + \cos(\gamma_r) \cos(\gamma_i) \cos(\phi_r - \phi_i)]. \quad (2) $$

In eq. (2), $\tau_0$ is the rupture duration, $L$ is the rupture length and $V_i$ is the velocity of the considered phase (P or S). The angles $\gamma_r$ and $\phi_r$ give the dip and azimuth of the rupture, while $\gamma_i$ and $\phi_i$ are the

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Variation in $\tau_i$ as a function of $\cos\theta_i/V_i$, following eq. (2) for the rupture orientation leading to the minimum misfit. Circles indicate measurements from P waves, while squares are measurements from S waves. The slope of the best fit (dashed grey lines) gives an estimate of $L$, while the intercept provides an estimate of $\tau_0$.</figcaption></figure>

take-off angle and azimuth of the ray to station $i$. Here, we aim to recover the rupture direction ($\gamma_r, \phi_r$) as well as its characteristics ($L, \tau_0$). To do this, we employ a grid search approach over all possible values of $\gamma_r$ and $\phi_r$ by a step angle of 1°. For each tested pair of angles, eq. (2) becomes

$$ \tau_i^{th}(\gamma_r, \phi_r) = \tau_0 + \frac{L}{V_i} \cos(\theta_i), \quad (3) $$

where $(\theta_i)$, the angle between the rupture and the take-off direction for the $i$th station, is obtained by identification from eq. (2). We then estimate the best $\tau_0$ and $L$ from a simple linear least-squares fit minimizing $S(\gamma_r, \phi_r) = \sum_i ||\tau_i^{th}(\gamma_r, \phi_r) - \tau_i||_2$ (Fig. 8). Take-off angles were obtained from our new locations within the 3-D model, and we use here astf recovered from both P and S waves such that $V_i$ can be the P- or S-wave velocities derived from the 3-D model at the source location. We infer the best rupture direction as the one given by the lowest value of $S(\gamma_r, \phi_r)$, which gives an azimuth of $N318^\circ E$ and dip of $\gamma_r = -36^\circ$, meaning that the rupture propagates towards the surface (Fig. S7, Supporting Information). From this rupture direction, the best-fitting estimates of $\tau_0$ and $L$ are 0.14 ± 0.004 s and 322 ± 45 m, respectively. Assuming $\tau_0 = L/V_r$, where $V_r$ is the rupture speed, we can deduce from the previous estimates that $V_r = 2301 \text{ m s}^{-1}$. We note that the azimuth, $N318^\circ E$, of the rupture is very close to the direction of the first nodal plane inferred from this earthquake, $N304^\circ E$ (Table 2). When we compute the dot product between the rupture direction and the normal vectors of the two nodal planes, we identify that the first nodal plane with a strike direction $N124^\circ E$ actually has a dot product three times smaller than the second nodal plane with a strike direction $N33^\circ E$. The directivity of the largest event in the sequence thus indicates that the fault plane corresponds to the plane in agreement with the main direction of relocated small earthquakes near the main shock and in a direction almost NW-SE (Fig. S8, Supporting Information).

As a consistency check, we can interpret the corner frequency as derived from the spectrum fitting in the previous section in terms of the source dimension. The relation between corner frequency and rupture length for unilateral rupture taken from (Savage 1972) gives
$L = 1.0 \times \beta/(\pi f_c)$, where we have considered that the ratio $V_R/\beta = 0.7$. We find $L=210$ m from $f_c = 5$ Hz for the largest event, which is, slightly lower, but still in agreement with the dimension deduced from the directivity analysis.

## 5 GROUND MOTION AND INTENSITY MODELLING

To link inverted source features at depth and inferred intensity maps at the surface, we carried out forward modelling of the ground motion based on the retrieved source parameters to see if the computed source properties and the derived 3-D velocity model are both in agreement with the observed intensities. To compute the seismic wavefield caused by the 2021 January earthquake, we employ the same finite-difference technique as in Section 4.1. Based on the similar corner frequency for the 2021 January and June events, we set the source duration at 0.14 s, and make use of the 3-D velocity model created in Section 2.1. The highest resolved frequency, according to the parameter settings, is 20 Hz. To assess the impact of the medium on the recovered waveforms, we also run a comparable simulation using the 1-D reference velocity model rather than the derived 3-D model. We convert the simulated velocity seismograms into seismic intensities using the maximum amplitude in the horizontal direction and the relation from Caprio et al. (2015) for each point of the grid that is located at the surface. We report the observed intensities for each map (Fig. 9), which were obtained from macroseismic data or computed from seismic instruments situated in the study area and converted to intensities using the same relation as for the simulated ones (Caprio et al. 2015).

We note that when we use the reference 1-D velocity model, the intensities we obtain are too low in comparison to the observed intensities. The largest computed intensities, which reflect the event's focal mechanism, are seen in four lobes surrounding the epicentre. The largest intensities in the 3-D model have a shape that is somewhat similar to the 1-D reference model, but the intensities are greater because the shallowest layers have less rigidity. As a result of the 3-D variations in the velocity model, we also observe some local effects. Notably, we can see some amplification along the Robertsau fault in the south. Overall, we see that the computed intensities in the 3-D case are higher and more consistent with the observed ones.

Our simulation is not intended to fully capture the complexity of the waveforms at high frequency. In fact, the simulated shaking can be altered by details of the source slip process, small size velocity perturbations and local site effects (which may be quantified by variation of the $V_s$30). However, we find that our method accurately reproduces the shaking intensity level without considering such refinements, and that the modelled intensities are quite accurate. This comparison shows how our source, location and velocity model results can all be interpreted within a unified framework.

## 6 DISCUSSION

### 6.1 Impact of the 3-D velocity structure

Our findings suggest that important 3-D velocity variations in the Rhine Graben must be taken into consideration, particularly for the location of events near the geothermal site. By applying an appropriate weighting of the various seismic wave phases based on their knowledge of the geological structural unit at each station site, the BCSF-Rnass operators reduce these effects. However, this

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Modelling of the seismic intensity produced by the 2021 January earthquake (background colours). Black lines show faults included in the GEORG database, and brown lines indicate simplified well trajectories. The epicentre location is depicted with a red star. Macroseismic intensities are displayed as coloured circles, and are estimated as average localities whose administrative contours are shown as thin white lines. The colour triangles refer to intensity derived from waveforms. Top: simulation performed in the reference 1-D velocity model. Bottom: simulation performed with the derived 3-D velocity model.</figcaption></figure>

method only uses some of all available phase picks for location and still produces some arrival times that are poorly predicted.

Here, we demonstrate that, for both *P* and *S* waves, all traveltime picks for all stations within regional distances, including those close to the epicentres in the sedimentary basin, can be reconciled using a 3-D velocity model at the regional scale. The 3-D model's effects go beyond just the earthquake location; they also have an impact on the computation of the moment magnitude and focal mechanism. In fact, the event's location can alter estimates of its depth and magnitude. Second, if the low velocity (low rigidity) of the subsurface material is not taken into account when estimating a moment magnitude based on the displacement recorded in the graben with a thick sedimentary cover, the magnitude may be overestimated. The evolution of the TLS, which is based on a magnitude threshold in this case, as well as our ability to understand the evolution of the
reservoir from the analysis of the seismicity, are both impacted by
these effects.

## 6.2 Identification of the active fault plane

Close to the wells, the orientation of the fault plane that slipped during the largest event is approximately N318E (N304E from the moment tensor inversion), in agreement with the orientation of the earthquake cloud obtained from precise double-difference re-location (Schmittbuhl et al. 2021) and Fig. S8 in the Supporting Information. We tested if this orientation is compatible with the stress field of the area. The orientation of the principal horizontal stress is deduced from the two closest well measurements located in Soultz-sous-Forêts and in Rittershoffen. Both sites are geothermal reservoirs located 30 km north of the GEOVEN reservoir.

According to Cornet et al. (2007), the maximum horizontal stress ($S_H$) at the Soultz-sous-Forêts geothermal site is oriented N170E. In the Rittershoffen reservoir, Azzola et al. (2019) and Hehn et al. (2016) found that the direction of $S_H$ changes with depth and that the granitic basement can reach N155E. Assuming that the principal stress directions at the location of the studied earthquakes are similar to these two cases, and hypothesizing a friction coefficient $\mu = 0.8$, we found that the optimally oriented fault planes are N325E and N015E for $S_H$ oriented N170E and N310E and N0E for $S_H$ oriented N155E. The fault plane identified here, N318E, has an orientation very close to the N325E and N310E optimal planes, indicating that it is the easiest to destabilize. Our results imply that identifying the active fault plane is best captured by the precise earthquake location and the accurate modelling of the source of the larger events, as done here. This task can only be achieved with a sufficiently dense network of equipment close to the earthquake locations.

## 6.3 Impact of the source parameters on the intensities

### 6.3.1 Earthquake magnitude

Our results show that even if the induced events analysed in this study have a moment magnitude $M_w$ only slightly larger than 3.0, their location, at the centre of a densely populated area, 4 km deep and in a sedimentary basin, still results in widespread shaking (up to intensity V) and some damage, as was also seen in a manner similar to that during the 2006 Basel earthquake (Deichmann & Giardini 2009). The location and source estimates of these four earthquakes help us to interpret the earthquake intensities and discuss their differences in pattern, notably for the three earthquakes that took place close to the wells because they share nearly identical locations and mechanisms. For these three events, we can test, at common sites, how the intensity changed between each shock. For all sites located less than 10 km from the epicentre, we computed the difference in intensity (as estimated from macroseismic data or derived from PGV) between the 2021 June earthquake and the other two earthquakes of 2020 December and 2021 January. The average intensity difference is −0.28 for the December earthquake and −0.44 for the 2021 January earthquake. We can check whether these differences in average intensity match the measured moment magnitude of these events. Indeed, we can convert the variation in intensity, $\Delta I$, into a difference in moment magnitude, $\Delta M_w$, using the ground motion prediction equation of Bakun & Scotti (2006) derived for the Rhine Graben, at a common site and at a similar epicentral distance, using $\Delta M_w = \Delta I/1.27$. This results in moment magnitude differences between the earthquakes in 2021 June and

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Variation in the intensity between the 2021 June earthquake and the 2021 January earthquake (from macroseismic data, circles and inferred from PGV, triangles). The average difference in intensity has been removed. The arrow shows the epicentre location of the June event and points in the inferred rupture direction.</figcaption></figure>

2020 December and between those in 2021 June and 2021 January, of −0.22 and −0.34, respectively, in accordance with the computed magnitudes (see Table 2).

### 6.3.2 Directivity

We also tested whether the rupture directivity of the 2021 June earthquake affected the recovered ground shaking associated with this event. This directivity effect manifests as an azimuthal modulation of the intensity amplitude. We computed the difference in intensity for each common site between the 2021 June event and the 2021 January earthquake (Fig. 10). Because the magnitude of the June event is larger than that of the January event, we expected the difference in intensity at all sites to be positive. To correct for this global effect, we removed the average intensity difference from each estimate to reveal only possible azimuthal variations. Because both tested earthquakes are located very close to each other and share almost the same focal mechanism, we expect that all the perceived differences in shaking at similar locations originate from source differences between the two events. We expect the difference in intensity at all sites to be positive since the magnitude of the June event was greater than the magnitude of the January event. We eliminated the average intensity difference from each estimate as a way to correct for this global effect and to only reveal potential azimuthal variation. We can interpret the observed intensity variations as the result of the rupture directivity of this larger event if we assume that the directivity effect of the 2021 January earthquake is small, which is reasonable given its smaller magnitude compared to the June 2021 earthquake. We note that almost all intensities show a consistent pattern of greater ground shaking to the NW of the epicentre (Fig. 10). Only lower amplitudes are seen in the opposite direction. The direction of rupture previously predicted for the earthquake in 2021 June corresponds to the direction where enhanced intensities are observed. This agreement between the direction of the rupture and the direction of the increased intensities demonstrates how the shaking caused by these events is being modified by such effects.
### 6.3.3 Stress drop and 3-D effects

The magnitude difference between the 2019 earthquake and the 2021 January event is relatively small, but the 2019 earthquake produced much lower intensities. For this observation, we can offer two main explanations. First, it is clear that the two events did not take place in the same place. This difference in location could highlight the importance of 3-D effects of the geological structures dependent on the source location, in modulating the severity of the perceived shakings. Indeed, the geographical variation in the thickness of the sediment cover in the area and the presence of faults that induce material contrasts (see Fig. S2, Supporting Information) can lead to local amplification or reduction in ground shaking depending on the direction of the incoming wavefield (Graves et al. 1998; Frankel et al. 2009). The 2019 earthquake's lower stress drop is another explanation for why it produced less ground motion than the 2020 earthquake. Our spectral analysis shows that the corner frequencies of both earthquakes are similar, suggesting that the rupture durations of both events are comparable. Since the first event's moment is slightly lower, its moment rate is also slightly lower over a comparable duration, making it a less energetic event with a smaller stress drop. As previously observed for injection-induced events, lower stress drop events have been linked to lower intensities (Hough 2014).

## 6.4 TLS design based on ground motion modelling

Our results illustrate how the perceived ground shaking related to injection-induced earthquakes in urban environments can change depending on the source and location of the earthquakes. It emphasizes the fact that the ground motion variability can be accounted for by geological factors and earthquake properties that can both have an impact on the perceived shaking. All these factors must be taken into account in the modelling of possible scenarios related to the development of a geothermal project (e.g. Ripperger et al. 2009). Identifying the possible faults and future earthquake locations, estimating the variability in earthquake source properties, including possible rupture directivity effects, and modelling the resulting wavefield for each scenario can then be a beneficial task. Indeed, such an approach would make it possible to estimate, even before the start of a geothermal project, the critical zones that might possibly suffer from the largest surface vibrations and thus to adapt the design of the project or even decide its continuation if the results of these simulations are unfavourable. For example, such a strategy has been successfully applied to investigate the effect of induced seismicity near Munich, Germany (Keil et al. 2022). In our study, because of the dense population all around the epicentres and the magnitude of the analysed events, the distribution of the population is not an essential factor to consider in these scenarios. However, depending on the project location, such information should also be used in establishing the outcome of such simulations, as illustrated by Schultz et al. (2021). Our results therefore encourage basing the adaptative-TLS system on ground shaking rather than on magnitudes. The combined effect of magnitude, focal mechanism and possibly rupture direction are all key aspects to consider as well as the development of realistic velocity models, especially for shallow unconsolidated sedimentary layers.

## 7 CONCLUSION

We have shown that, taking into account a realistic velocity model for the region, we are able to recover the source parameters of the

four most important earthquakes that occurred during the Strasbourg induced seismicity sequence. Modelling the ground motion from these parameters shows that the predicted intensity is in agreement with that observed in the vicinity of the geothermal site. As the population's perception of the severity of shaking may impact the development of geothermal projects, we propose that scenarios that consider the modelling of the wavefield produced by induced earthquakes be included in TLS.

## ACKNOWLEDGMENTS

We thank Zacharie Duputel and Luis Rivera for discussion on the analysis of the directivity. We thank Andrea Bruestle from LGRD for German's macroseismic data compilation. We thank the two anonymous reviewers and the editors that provided numerous comments that significantly improved the discussion between the source process and the observed intensities. This work was carried out in the framework of the Interdisciplinary Thematic Institute GeoT as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm. It was supported by ANR PrESENCE (ANR-21-CE05-0033-01), IdEx Unistra (ANR-10-IDEX-0002), and by SFRI-STRAT'US project (ANR ANR-20-SFRI-001) under the framework of the French Investments for the Future Program.

## DATA AVAILABILITY

All waveforms used in this study are available through the RESIF (RESIF 1995) and RASPBERRY-SHAKE (Raspberry Shake 2016) interfaces. The earthquake catalogue of the BCSF-Rénass, waveforms and the 3-D velocity model built in this study are all available in a dedicated repository at https://cdgp.u-strasbg.fr. Some of the data processing is accomplished using Obspy software (Krischer et al. 2015) and several figures have been made with GMT (Wessel et al. 2013). Some figures also use map tiles by Stamen Design under CC BY 3.0 and data by OpenStreetMap under ODbL.

## REFERENCES

* Abercrombie, R.E., Poli, P. & Bannister, S., 2017. earthquake directivity, orientation, and stress drop within the subducting plate at the hikurangi margin, new zealand, *J. geophys. Res.: Solid Earth*, **122**(12), 10–176.

* Aki, K. & Richards, P.G., 2002. *Quantitative Seismology*, University Science Books, 2nd edn.

* Azzola, J., Valley, B., Schmittbuhl, J. & Genter, A., 2019. Stress characterization and temporal evolution of borehole failure at the rittershoffen geothermal project, *Solid Earth*, **10**(4), 1155–1180.

* Bakun, W.H. & Scotti, O., 2006. Regional intensity attenuation models for france and the estimation of magnitude and location of historical earthquakes, *Geophys. J. Int.*, **164**(3), 596–610.

* Bertero, M., Bindi, D., Boccacci, P., Cattaneo, M., Eva, C. & Lanza, V., 1997. Application of the projected landweber method to the estimation of the source time function in seismology, *Inverse Prob.*, **13**(2), 465, doi:10.1088/0266-5611/13/2/017.

* Bohnhoff, M., Malin, P., ter Heege, J., Deflandre, J.-P. & Sicking, C., 2018. Suggested best practice for seismic monitoring and characterization of non-conventional reservoirs, *First Break*, **36**(2), 59–64.

* Bommer, J.J. & Crowley, H., 2006. The influence of ground-motion variability in earthquake loss modelling, *Bull. Earthq. Eng.*, **4**(3), 231–248.

* Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, *J. geophys. Res.*, **75**(26), 4997–5009.
Caprio, M., Tarigan, B., Worden, C.B., Wiemer, S. & Wald, D.J., 2015. Ground motion to intensity conversion equations (gmices): A global relationship and evaluation of regional dependency, *Bull. seism. Soc. Am.*, **105**(3), 1476–1490.

Charléty, J., Cuenot, N., Dorbath, C. & Dorbath, L., 2006. Tomographic study of the seismic velocity at the soultz-sous-forêts eggs/hdr site, *Geothermics*, **35**(5–6), 532–543.

Cornet, F. & Jianmin, Y., 1995. Analysis of induced seismicity for stress field determination and pore pressure mapping, in *Mechanics Problems in Geodynamics Part I*, pp. 677–700, Springer.

Cornet, F., Helm, J., Poitrenaud, H. & Etchecopar, A., 1997. Seismic and aseismic slips induced by large-scale fluid injections, in *Seismicity Associated with Mines, Reservoirs and Fluid Injections*, pp. 563–583, Springer.

Cornet, F., Bérard, T. & Bourouis, S., 2007. How close to failure is a granite rock mass at a 5 km depth?, *Int. J. Rock Mech. Mining Sci.*, **44**(1), 47–66.

Cornet, F.H., 2016. Seismic and aseismic motions generated by fluid injections, *Geomech. Energy Environ.*, **5**, 42–54.

Deichmann, N. & Giardini, D., 2009. Earthquakes induced by the stimulation of an enhanced geothermal system below basel (switzerland), *Seismol. Res. Lett.*, **80**(5), 784–798.

Foulger, G.R., Wilson, M.P., Gluyas, J.G., Julian, B.R. & Davies, R.J., 2018. Global review of human-induced earthquakes, *Earth-Sci. Rev.*, **178**, 438–514.

Frankel, A., Stephenson, W. & Carver, D., 2009. Sedimentary basin effects in Seattle, Washington: ground-motion observations and 3d simulations, *Bull. seism. Soc. Am.*, **99**(3), 1579–1611.

Freymark, J., et al., 2020. 3d-urg: 3d gravity constrained structural model of the upper rhine graben, *GFZ Data Services*, doi:10.5880/GFZ.4.5.20.20.004.

Graves, R.W., Pitarka, A. & Somerville, P.G., 1998. Ground-motion amplification in the santa monica area: Effects of shallow basin-edge structure, *Bull. seism. Soc. Am.*, **88**(5), 1224–1242.

Grigoli, F. et al., 2017. Current challenges in monitoring, discrimination, and management of induced seismicity related to underground industrial activities: A european perspective, *Rev. Geophys.*, **55**(2), 310–340.

Grünthal, G., 1998. *European macroseismic scale 1998*, Tech. Rep., European Seismological Commission (ESC).

Guglielmi, Y., Cappa, F., Avouac, J.-P., Henry, P. & Elsworth, D., 2015. Seismicity triggered by fluid injection–induced aseismic slip, *Science*, **348**(6240), 1224–1226.

Hanks, T.C. & Kanamori, H., 1979. A moment magnitude scale, *J. geophys. Res.: Solid Earth*, **84**(B5), 2348–2350.

Haskell, N., 1964. Total energy and energy spectral density of elastic wave radiation from propagating faults, *Bull. seism. Soc. Am.*, **54**(6A), 1811–1841.

Hehn, R., Genter, A., Vidal, J. & Baujard, C., 2016. Stress field rotation in the EGS well GRT-1 (Rittershoffen, France), in *Proceedings of European Geothermal Congress*, European Geothermal Energy Council, ISBN: 978-2-9601946-0-9.

Hopp, C., Sewell, S., Mroczek, S., Savage, M. & Townend, J., 2019. Seismic response to injection well stimulation in a high-temperature, high-permeability reservoir, *Geochem. Geophys. Geosyst.*, **20**(6), 2848–2871.

Hough, S.E., 2014. Shaking from injection-induced earthquakes in the central and eastern united states, *Bull. seism. Soc. Am.*, **104**(5), 2619–2626.

Keil, S., Wassermann, J. & Megies, T., 2022. Estimation of ground motion due to induced seismicity at a geothermal power plant near munich, Germany, using numerical simulations, *Geothermics*, **106**, 102577, doi: 10.1016/j.geothermics.2022.102577.

Krischer, L., Megies, T., Barsch, R., Beyreuther, M., Lecocq, T., Caudron, C. & Wassermann, J., 2015. Obspy: a bridge for seismology into the scientific python ecosystem, *Comput. Sci. Discov.*, **8**(1), 014003, doi:10.1088/1749-4699/8/1/014003.

Kwiatek, G. et al., 2019. Controlling fluid-induced seismicity during a 6.1-km-deep geothermal stimulation in Finland, *Sci. Adv.*, **5**(5), eaav7224, doi: 10.1126/sciadv.aav7224.

Lanza, V., Spallarossa, D., Cattaneo, M., Bindi, D. & Augliera, P., 1999. Source parameters of small events using constrained deconvolution with empirical green's functions, *Geophys. J. Int.*, **137**(3), 651–662.

Lengliné, O., Boubacar, M. & Schmittbuhl, J., 2017. Seismicity related to the hydraulic stimulation of grt1, rittershoffen, france, *Geophys. J. Int.*, **208**(3), 1704–1715.

Lomax, A., Zollo, A., Capuano, P. & Virieux, J., 2001. Precise, absolute earthquake location under somma-vesuvius volcano using a new three-dimensional velocity model, *Geophys. J. Int.*, **146**(2), 313–331.

Maeda, T., Takemura, S. & Furumura, T., 2017. Openswpc: an open-source integrated parallel simulation code for modeling seismic wave propagation in 3d heterogeneous viscoelastic media, *Earth Planets Space*, **69**(1), 1–20.

Maurer, V., Gaucher, E., Grunberg, M., Koepke, R., Pestourie, R. & Cuenot, N., 2020. Seismicity induced during the development of the rittershoffen geothermal field, france, *Geothermal Energy*, **8**(1), 1–31.

McClure, M.W., 2012. *Modeling and Characterization of Hydraulic Stimulation and Induced Seismicity in Geothermal and Shale Gas Reservoirs*, Stanford University.

Mignan, A., Broccardo, M., Wiemer, S. & Giardini, D., 2017. Induced seismicity closed-form traffic light system for actuarial decision-making during deep fluid injections, *Sci. Rep.*, **7**(1), 1–10.

Park, S. & Ishii, M., 2015. Inversion for rupture properties based upon 3-d directivity effect and application to deep earthquakes in the sea of okhotsk region, *Geophys. J. Int.*, **203**(2), 1011–1025.

Raspberry Shake, S.A., 2016. *Raspberry shake*, doi:10.7914/SN/AM.

RESIF, 1995. RESIF-RLBP french broad-band network, RESIF-RAP strong motion network and other seismic stations in metropolitan France, doi:10.17616/R37Q06.

Ripperger, J., Kästli, P., Fäh, D. & Giardini, D., 2009. Ground motion and macroseismic intensities of a seismic event related to geothermal reservoir stimulation below the city of basel-observations and modelling, *Geophys. J. Int.*, **179**(3), 1757–1771.

Rothé, J.-P. & Peterschmitt, E., 1950. *Etude séismique des explosions d'Haslach*, impr. Alsacienne, *Ann. Inst. Phys. Globe Strasbourg*, **5**(3), 3–28.

Satriano, C., 2021. *Sourcespec – earthquake source parameters from s-wave displacement spectra*, doi:10.5281/zenodo.6954238.

Savage, J., 1972. Relation of corner frequency to fault dimensions, *J. geophys. Res.*, **77**(20), 3788–3795.

Schmittbuhl, J., Lengliné, O., Cornet, F., Cuenot, N. & Genter, A., 2014. Induced seismicity in eggs reservoir: the creep route, *Geothermal Energy*, **2**(1), 1–13.

Schmittbuhl, J. et al., 2021. Induced and triggered seismicity below the city of strasbourg, france from november 2019 to january 2021, *C. R. Géosci.*, **353**(S1), 1–24.

Schoenball, M., Dorbath, L., Gaucher, E., Wellmann, J.F. & Kohl, T., 2014. Change of stress regime during geothermal reservoir stimulation, *Geophys. Res. Lett.*, **41**(4), 1163–1170.

Schultz, R., Beroza, G.C. & Ellsworth, W.L., 2021. A strategy for choosing red-light thresholds to manage hydraulic fracturing induced seismicity in North America, *J. geophys. Res.: Solid Earth*, **126**(12), e2021JB022340, doi:10.1029/2021JB022340.

Shapiro, S.A., Rothert, E., Rath, V. & Rindschwentner, J., 2002. Characterization of fluid transport properties of reservoirs using induced microseismicity, *Geophysics*, **67**(1), 212–220.

Sira, C. & Grunberg, M., 2021. Rapport macrosismique - séisme induit de la Wantzenau (Bas-Rhin), 26 juin 2021 à 5h00 locale, *Magnitude 3.9 ML(RENASS)*, Intensité maximale V (EMS98), Tech. rep., BCSF-Rénass, Strasbourg, France.

Sira, C., Cara, M., Schlupp, A., Masson, F., Schaming, M. & Mendel, V., 2021a. 1921-2021: 100 years of macroseismic studies at bcsf, *C. R. Géosci.*, **353**(S1), 23–51.

Sira, C., Mendel, V., Schaming, M., Cara, M., Masson, F. & Schlupp, A., 2021b. Collection des différents formulaires métropolitains d'enquête macrosismique des archives du bureau central sismologique français 1910 à décembre 2020., Tech. rep., BCSF-Rénass, Strasbourg, France.

Vallée, M., 2004. Stabilizing the empirical green function analysis: Development of the projected landweber method, *Bull. seism. Soc. Am.*, **94**(2), 394–409.
Wang, X. & Zhan, Z., 2020. Moving from 1-d to 3-d velocity model: automated waveform-based earthquake moment tensor inversion in the los angeles region, *Geophys. J. Int.*, **220**(1), 218-234.

Wei, S. et al., 2015. The 2012 brawley swarm triggered by injection-induced aseismic slip, *Earth planet. Sci. Lett.*, **422**, 115-125.

Wessel, P., Smith, W.H., Scharroo, R., Luis, J. & Wobbe, F., 2013. Generic mapping tools: improved version released, *Eos, Trans. Am. geophys. Un.*, **94**(45), 409-410.

Zang, A., Oye, V., Jousset, P., Deichmann, N., Gritto, R., McGarr, A., Majer, E. & Bruhn, D., 2014. Analysis of induced seismicity in geothermal reservoirs-an overview, *Geothermics*, **52**, 6-21.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Difference of *P*-wave velocity between the reference 1-D model and the constructed 3-D model at 2 km depth. The green triangle is the location of the GEOVEN wellheads and the blue triangles are the location of seismic stations used in this study.

**Figure S2.** Traveltime residuals for *P* wave (left) and *S* wave (right) computed for the 1-D reference model (blue), and for the derived regional 3-D model (red). The large delays on *S*-wave traveltimes visible for the 1-D model mostly are much reduced using the 3-D model.

**Figure S3.** Histogram showing the distribution of the depth of scatter points drawn from the pdf of the location of each event. It highlights the similar resolved depth for all events. For comparison, the coloured dashed lines show the depth of the same event in the BCSF-Rénass catalogue. Colours refer to each earthquake similarly as in Fig. 2.

**Figure S4.** Intensity maps for the four studied events. For the 2019 event, the intensity map is from citizen testimonies received by the BCSF-Rénass and LGRB, averaged by localities. For the three others, intensities are estimated using communal and individuals forms for the French part and only individuals forms from localities in Germany. The red stars give the location of the epicentres.

**Figure S5.** Waveforms recorded at two sites (left: site in Reichstett, station R60B1, approximately 3 km from the reservoir, right: Strasbourg seismology museum, station STMU, approximately 9 km from the reservoir). All records start at the origin time of each event and last 15 s. The records at R60B1 are velocity in m s⁻¹, while at STMU are ground motion acceleration in m s⁻².

**Figure S6.** Example of the deconvolution process at station R461D. Left: the top panel shows the *P*-wave record of the EGF at station R461D. The middle panel shows the recording of the *P* wave as well at the same station for the 2021 June event, that is, the main event. On the bottom panel, the red curve shows again the main event waveform while the black line shows the recovered signal by convolution of the EGF signal with the apparent source time function displayed on the right. Right: evolution of the misfit function ϵ as a function of the total duration allowed for the deconvolution process, *T*. Here we observe a rapid decrease of ϵ up to *T* ≃ 0.1 s and after an almost constant level of ϵ. The dashed red curve indicates the fixed value of *T* based on this curve.

**Figure S7.** Variation of the misfit as a function of the orientation of the rupture direction. The misfit is the difference of the sum of square between the observed and computed τᵢ obtained for each orientation of the rupture direction. Negative values indicate an upward rupture. The lowest misfit is obtained for φᵣ = 318° and γᵣ = -36°.

**Figure S8.** Map of relocated events up to the 2021 June 28 (grey dots). Initial, absolute locations of the relocation process are taken from the BCSF-Rénass locations. Red stars are the relocations of the three largest event in the area. The black curve shows the simplified trajectory of the GT2 well. The details of the relocation process follows the same procedure as in Schmittbuhl *et al.* (2021).

**Table S1.** Quality check estimators of the moment tensor inversion. We report the shift in origin time, Δ*T*₀ for each event that best match the waveform and the average correlation coefficient, *CC* between synthetics and observed waveforms.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.