Rupture Process of the Mw 3.3 Earthquake in the St. Gallen 2013 Geothermal Reservoir, Switzerland

We analyze slip distribution and rupture kinematics of a Mw3.3 induced event that occurred in the St. Gallen geothermal reservoir (NE Switzerland) in 2013. We carry out a two‐step procedure: (1) path effects are deconvolved from the seismograms using an empirical Green's function, resulting in relative source time functions at all seismic stations; (2) the relative source time functions are back‐projected to the corresponding isochrones on the fault plane. Results reveal that the mainshock rupture propagates toward NNE from the hypocenter with an average velocity of 2,000 m/s. Spatiotemporal organization of foreshocks and aftershocks shows that the mainshock broke a previously less active portion of the fault and suggests that the aftershock sequence could be mainly driven by stress transfer. Applying this method in an operational environment could enable fast retrieval of seismic slip, allowing assessment of fault asperities and structures involved in the reservoir creation process.

In most models and analyses, earthquakes are considered either point sources or "penny-shaped" surfaces with homogeneous source properties. While these assumptions may be valid, details of individual  Heuberger et al., 2016). Focal mechanisms (Diehl et al., 2017) of the mainshock (black), and the M L 2.1 event (EGF1 in red) are represented. A smaller M L 1.4 event (EGF2 in orange) is also shown; its focal mechanism is not available: here we attribute the same mechanism as the EGF1 based on waveform similarity. Inset shows the location of the region. earthquakes are more complex. Due to the insufficient resolution of most monitoring networks, micro-earthquakes (M < 2) must be still approximated as dislocation on penny-shaped surfaces. But in case of larger events (M > 3) recorded on well-distributed local seismic stations, detailed slip patterns within the rupture area can be imaged. Previous studies of small to moderate earthquakes of natural origin showed that, similarly to large earthquakes, they exhibit rupture directivity (e.g., Abercrombie et al., 2017;Boatwright, 2007;McGuire, 2002) and spatial variability in slip distribution (e.g., Dreger et al., 2007;Ide, 2001;Kim et al., 2016;McGuire, 2004). Imaging and modeling the rupture kinematics of earthquakes can contribute to the understanding of the structures associated with rupture nucleation, propagation, slipped asperities, and seismic radiation (e.g., Mai, 2005). Moreover, the comparison of the kinematic slip image with the space-time distribution of earthquakes that occurred before and after can provide insights into the load processes leading to seismic failure on those structures (e.g., Stein, 2003). Slip distributions can also serve as an input or independent information for modeling. For instance, the slip distribution of larger events is required for models incorporating static stress change triggered by earthquakes, or hydromechanical models that require a relationship between slip and hydraulic aperture (e.g., Rutqvist, 2011).
In this study, we examine the slip distribution and the rupture kinematics of a M w 3.3 (M L 3.5) induced event that occurred in the St. Gallen geothermal reservoir (NE Switzerland, Figure 1) on 20 July 2013, which was widely felt in the region (Edwards et al., 2015). The industrial project was inspired by hydrothermal systems where hot water is extracted from a deep natural reservoir (e.g., Evans et al., 2012). To test and improve the connectivity between the well and the naturally fractured aquifers, for example, fault zones, a hydraulic test and two acid stimulations were carried out in July 2013, immediately followed by seismic events (Diehl et al., 2017;Wolfgramm et al., 2015). On 19 July, methane gas was detected in the borehole. To reduce the pressure buildup, drilling mud was pumped into the well, increasing significantly the micro-earthquake  Diehl et al., 2017), followed by hundreds of events in the next few months. The St. Gallen sequence ceased completely on 25 October 2013, following a production test (e.g., Moeck et al., 2015). Seismic monitoring was carried out with a local seismic network consisting of one short period borehole sensor at about 200-m depth (SGT00), five broadband surface stations (SGT01-SGT05), and additional ten short-period surface stations (SGT06-15; Kraft et al., 2013), from which SGT00 and SGT02-06 are selected for this study based on signal quality (Figures 1 and 2). Between July and November 2013, 347 earthquakes were located and subsequently finely relocated by Diehl et al. (2017). Most of the seismicity clusters in an area of roughly 0.5 × 0.4 km 2 , mainly coinciding with the supposed rupture area of the largest event, inferred from scaling laws (Edwards et al., 2015). The details of the rupture process of the main event are, however, poorly known. Our goal, here, is to give insight into the rupture kinematics of the mainshock (MSH; e.g., rupture directivity and velocity) and into the spatiotemporal organization of the MSH slip on the fault, relative to the seismic events that occurred before and after.

Method
To image the slip pattern of the main St. Gallen earthquake, we use a two-step procedure (Adinolfi et al., 2015;Stabile et al., 2012), which treats separately the problems of modeling the wave propagation and constructing the kinematic slip model. In the first step, the propagation effects are deconvolved from the recorded seismograms by applying an empirical Green's function (EGF), that is, a smaller event located close to the MSH and sharing a similar focal mechanism (e.g., Capuano et al., 1994;Courboulex et al., 1996;Hartzell, 1978;Irikura, 1983). The deconvolution process yields an apparent moment rate function at each station-carrying information on source directivity and moment release history-which we call "relative source time function" (RSTF). In the second step, the RSTFs are back-projected to the fault plane through an "isochrone back-projection" approach (IBP; Festa & Zollo, 2006). These two steps, along with the results, are detailed in the following two sections.

RSTFs and Rupture Directivity
Recordings for the MSH, and of the two selected EGFs (M L 2.1-EGF1-and M L 1.4-EGF2-see map in Figure 1) are shown in Figure 2. Visual comparison of the traces evidences that MSH and EGF waveforms have similar features, though MSH recordings are characterized by lower frequencies due to the extended source process. This also confirms the similar focal mechanism between MSH and EGF1 found by Diehl et al. (2017). No focal mechanism is available for EGF2, though waveform similarity indicates that it is comparable to EGF1. MSH signals present a small additional onset before the first arriving P wave (Figure 2, insets), which we interpret as the P wave from a small earthquake, occurring a few hundredths of seconds before the MSH.
Deconvolution of the EGF from the MSH cancels the common path and radiation pattern effects, as well as the instrumental response. The resulting RSTF represents the scalar moment release of the earthquake, distorted by the directivity effect due to the rupture propagation along the fault. Signal deconvolution is carried through the projected Landweber method (Vallée, 2004), a constrained iterative least-squares minimization approach which enforces nonnegativeness, causality, finite duration, conservation of moment ratio of the target event, and EGF (supporting information Figure S1). Traces are filtered between 0.1 and 15 Hz, that is, below the theoretical corner frequency of the two EGFs (16 and 27 Hz, respectively, assuming a stress drop of 3 MPa and a shear wave velocity of 3 km/s; Edwards et al., 2015). The deconvolution is performed separately for the P waves (using vertical components) and S waves (using north-south and east-west components independently). In the following, we shall call the resulting P wave RSTFs "P RSTFs" and the S wave RSTFs with north-south and east-west components "Sn RSTFs," and "Se RSTFs," respectively. Before deconvolution, MSH and EGF traces need to be aligned. This cannot be done independently for each station, since it could produce incoherent results. We follow instead a master event approach (e.g, Grigoli et al., 2016) and search for the best centroid of MSH moment release, with respect to a fixed EGF location, as the point on the fault that better fits the relative shift between MSH and EGF P wave traces ( Figures S2 and S3). Centroid alignment avoids significant distortions in the RSTFs shape. Finally, to pinpoint the rupture start, the start of the EGFs is handpicked. Figure S4 shows RSTFs obtained using EGF1 and EGF2: the overall similarity between P, Sn, and Se RSTF for both EGFs validates the stability of the deconvolution process and the reliability of the retrieved RSTFs. The early onset on the MSH signal (shown in Figure 2) appears, after deconvolution, on the RSTFs for both EGF1 and EGF2. To image exclusively the main moment release, we pick the start of the highest moment release (vertical bars in Figure S4) and use only the remaining part of the RSTFs ( Figure S5) to construct the kinematic slip model. Figure 3 shows the final RSTFs obtained from deconvolution with EGF1, as a function of station azimuth from the rupture strike (N32), compared to the theoretical apparent source duration for a horizontal (i.e., along-strike) rupture propagation on a vertical fault (Haskell, 1964): where L is the rupture length (200 m), Vr is the rupture speed (1,500, 2,000, or 2,500 m/s), and c is the P or S wave velocity (5 and 3.5 km/s, respectively). This simple analysis indicates a NNE directivity of the rupture wave RSTFs, obtained from deconvolution of NS ("Sn RSTF") and EW ("Se RSTF") components, respectively. Curve colors correspond to station colors in Figure 1. Sn RSTF at station SGT04 is not used due to its low quality. Gray curves indicate the theoretical apparent source duration (using the standard directivity formula; Haskell, 1964) for a rupture length of 200 m in N32 direction and a rupture speed (Vr) of 1,500, 2,000, or 2,500 m/s; P and S wave velocities are fixed to 5 and 3.5 km/s, respectively.
with a speed between 1,500 and 2,000 m/s. The impossibility of finding a unique rupture speed implies that the rupture process is more complex than a linear propagation, as we shall see in the following section.

Kinematic Slip Model From IBP
In the second step we seek to back-project the RSTFs on the fault plane and reconstruct, through an averaging procedure, the space-time distribution of fault slip. We set up a 600-m × 600-m subvertical fault (dip of 84 • ) with 5-m × 5-m uniform grid cells, centered around the hypocenter of the MSH (47.421073 • , 9.319911 • , depth 4.332 km) (Diehl et al., 2017). Using the P and S wave velocity model from Diehl et al. (2017), and fixing the same constant rupture velocity in both along-dip and along-strike directions (i.e., circular rupture), the P, Sn, and Se RSTFs are back-projected over the rupture isochrones, that is, points on the fault plane that are "seen" as rupturing at the same time from the relative station position (Bernard & Madariaga, 1984;Spudich & Frazer, 1984). Without any a priori information on the moment rate on the fault, the RSTFs are uniformly distributed (back-projected) on the corresponding isochrones, and the resulting moment rate map (IBP image) is the spatial average on the fault of the back-projected RSTFs. The quality of the obtained solution is computed as the mean RSTF residual, using an L1 norm for the difference between the original and the expected RSTFs, reconstructed by solving the forward problem. As discussed in Festa and Zollo (2006), the obtained IBP map is a defocused image of the actual moment rate distribution on the fault, due to the convolutional effect introduced by the isochrone geometry. An imperfect-yet effective-deconvolution is achieved through the "restarting" iterative procedure proposed by Festa and Zollo (2006): A new back-projection is performed using the previous mean IBP image as spatial weighting. With this approach, moment rate is no longer spread uniformly on the isochrones but concentrated on the zones indicated by the previous mean IBP solution ( Figures S7-S9), resulting into a more focused IBP map and lower RSTFs residuals. The process is iterated while the mean RSTF residual decreases significantly (see the supporting information for more details and testing of the procedure for a synthetic case). We apply a Gaussian smoothing (with standard deviation of 0.5 grid cells) on every iteration to avoid too rough intermediate and final solutions, while keeping the same general properties of the moment rate pattern. Finally, moment rate is converted to seismic moment by assuming uniform rise time on the entire fault (equation (S22)) and to fault slip using Aki's seismic moment formula (equation (S24)).
Joint back-projection of P and S RSTFs, repeated for different rupture velocities, provides an optimal rupture velocity of 2,000 m/s ( Figure S13) and reveals a nonuniform slip pattern with significant rupture directivity ( Figure 4). The rupture initially propagated concentrically with about 10 cm of slip around the hypocenter, then moved toward the NNE activating an asperity about 150 m away from the hypocenter as well as a fairly separated patch below the main slip asperity with about 5 cm of slip. The bulk of the slip pattern is more compact than a circular rupture area corresponding to the stress drop values estimated by Edwards et al. (2015), which would suggest a higher local stress drop value (15-30 MPa, dark and light gray dotted circles, respectively, in Figure 4).

Discussion
According to the catalog produced by Diehl et al. (2017; magnitude of completeness: M L 0.0), the main event in the St. Gallen reservoir was preceded by 61 earthquakes (−0.9 < M L < 2.1) and was followed by 282 events (−1.2 < M L < 1.7). For simplicity, we call these events "foreshocks" and "aftershocks", respectively, referring solely to their origin time relative to the MSH. The space-time seismicity pattern is shown in Figure 5, with symbol sizes corresponding to rupture area, assuming a stress drop of 3 MPa and using Eshelby's formula (equation (S27)); T0 is the timing of the MSH. Early events (6 days to ∼ 10 hr before T0, Figure 5a) occur to the SSW of the MSH rupture. Their locations might be explained by stress concentration on the edges of a locked part of the fault, as suggested by Diehl et al. (2017). Note that there is no direct contact between the seismicity cloud and the openhole section of the well, but hydraulic connection must be already established as seismicity reacts promptly to stimulation procedures (Diehl et al., 2017). Later foreshocks (∼ 10 hr to 3 min before T0) migrate toward the middle of the MSH slip patch, in what could be interpreted as a slow nucleation phase: The main rupture initiates at the location of the largest foreshock (M L 2.1, EGF1 in this study), and the bulk of slip coincides with the foreshock location. On the other hand, the slip to the NNE appears to be complementary to the foreshock distribution. Early aftershocks (up to 5 days after T0) are mainly located around the MSH slip area (Figure 5b). Later aftershocks correspond to additional activities in the well (e.g., "fishing" for lost objects, cleaning, and production) and expand over a larger fault area, from SSW to NNE, with a significantly higher activity toward the NNE (Figure 5c). Diehl et al. (2017) already hypothesized that aftershock migration toward the NNE might have happened due to the gas kick that could alter hydraulic conditions or due to the main event that caused changes in hydraulic diffusivity by breaking a previously impermeable seal (possibly a NS striking fault). Our results suggest the plausibility of the latter hypothesis: it seems probable that the updip and NNE propagating MSH released a previously less active segment of the fault where significant seismicity occurred after the main event.
The observation of complementary organization of MSH slip and aftershocks is not new in the context of natural events, for example, M L 5.8 Coyote Lake (1979), California; M L 6.2 Morgan Hill (1984), California (Bakun et al., 1986); M s 7.3 Borah Peak (1983), Idaho; M s 6.0 North Palm Springs (1986), California (Mendoza & Hartzell, 1988); M s 8.1 Michoacan (1985), Mexico (Das, 2003;Mendoza & Hartzell, 1989); M w 7.3 Landers (1992), California (Das, 2003); M w 8.1 Antarctic Plate (1998;Das, 2003); M w 6.3 L'Aquila (2009), Italy (Chiarabba et al., 2009); M w 2.9 Irpinia (2008), Italy (Stabile et al., 2012); and M w 9.0 Tohoku-Oki (2011), Japan (Kato & Igarashi, 2012). Complementary patterns between early aftershocks and coseismic slip have also been recently observed for the induced M w 4.0 earthquake near Guthrie, Oklahoma (Wu et al., 2019). These observations are often explained by the idea that aftershocks occur where stress level is high after the MSH. These areas tend to be near the edges of the coseismic slip due to the increased loading on unbroken or partially broken barriers or in the region of rapid transition from high to low slip (Das, 2003;Mendoza & Hartzell, 1988). In the case of fluid-induced earthquakes, evolution of pore pressure due to fluid injection must play a role in the spatiotemporal occurrence of aftershocks, in addition to static or dynamic stress transfer from the MSH. To understand the role of the governing physical processes in the St. Gallen sequence, detailed geomechanical modeling is being carried out. It is worth stressing here the importance of a detailed MSH slip model: Using a homogeneous circular or rectangular rupture area might lead to misinterpretation of reservoir size and localized seismic slip ( Figure S15).
The inferred rupture model indicates directivity toward the injection well (Figures 1 and 5c), that is, "backward" rupture, following the terminology of Lui and Huang (2019). According to the 1-D numerical models of Dempsey and Suckale (2016), this would imply a high pressurization/criticality ratio (i.e., the ratio between injection pressure and initial fault stress) for the St. Gallen geothermal site.

Conclusions
Small earthquakes (M < 4) reveal their rupture complexity when observed close and when path effects are properly accounted for (e.g., by EGFs).
By retrieving RSTFs and back-projecting them to fault isochrones, we unraveled a two-stage rupture process for the M w 3.3 induced St. Gallen earthquake: (1) The bulk of the MSH slip (∼ 10 cm) takes place around the hypocenter, suggesting higher local stress drop (15-30 MPa), and (2) the rupture further propagates ∼ 150 m toward NNW at 2,000 m/s. Spatiotemporal organization of microseismicity before and after the MSH shows that (1) foreshocks migrated into the bulk of MSH slip area a few hours before the main rupture, evidencing a possible slow nucleation process; (2) the MSH broke-and possibly unsealed-a previously less active portion of the fault; and (3) aftershocks distribute around MSH slip area, suggesting that stress transfer from the MSH played a major role in the sequence.
The method presented here requires a high-quality-yet standard-microseismic monitoring system and low computational power, making it possible to implement the algorithm in an operational environment. This would allow for rapid assessment of fault asperities and structures involved in the reservoir creation process. Detailed slip models for fluid-induced events could further help in understanding the relationship between the location of seismic slip and fracture connectivity, retrieved through tracer tests (Hirschberg et al., 2015). This, in turn, might help in the future to define fluid paths in the reservoir prior to tracer tests and might give some hints about the optimal location of new boreholes. information. Conversion of latitude-longitude to Swiss coordinates was done using http://tool-online. com/en/coordinate-converter.php. The map of the study area was compiled using the following resources: map tiles by Stamen Design (https:// stamen.com), under CC BY 3.0., settlements by OpenStreetMap (https://www.openstreetmap.org), under the Open Database License, and roads and railways by https:// opendata.swiss, under Open use BY ASK. Data analysis was implemented using Python and the Numpy (Oliphant, 2007) and ObsPy (Krischer et al., 2015) libraries. Figures were prepared using Matplotlib (Hunter, 2007). We are grateful for discussions with Tobias Diehl, Toni Kraft, Martin Mai, Antonella Orefice, Tony Alfredo Stabile, and Martin Vallée during the preparation of the manuscript. We are also grateful to two anonymous reviewers and to Editor Gavin Hayes for their comments and suggestions which helped improve the manuscript.