Nonsystematic Rupture Directivity of Geothermal Energy Induced Microseismicity in Helsinki, Finland

The rupture behavior of microseismicity in fluid‐injection settings with low fault stresses is generally believed to be controlled by the pore pressure, including a tendency of the larger induced earthquakes to rupture into the perturbed volume toward the injection well, implying a degree of predictability. Here, we examine directivity patterns to identify fault planes and rupture directions of the 21 largest earthquakes (local magnitudes, ML ${M}_{L}$ 1.3–1.9) recorded during the 2018 St1 Deep Heat geothermal project near Helsinki, Finland. We use the Empirical Green's Function technique to retrieve per‐seismic‐station corner frequencies, earthquake durations, and directivity trends. After combining the directivity trends with focal mechanisms calculated using principle component analysis, we resolve rupture planes and rupture directions of 10 events. In contrast to studies of induced events at other sites, we find that one event rupture toward, two rupture away, and the remaining rupture parallel to the well. Furthermore, we find that the events prefer mode II failures rather than mode III. These observations provide new constraints for mechanical models of rupture growth in pore‐pressure dominated settings.


Data
In 2018, the 6.1 km deep OTN-3 geothermal well in Espoo, Finland, was hydraulically stimulated over 49 days (Hillers et al., 2020;Kwiatek et al., 2019). Two geophone arrays were deployed during the operation (see Figure 1). The nearby OTN-2 well contained ten high-frequency geophones optimal for microseismic detection (e.g., Kwiatek et al., 2019), while the shallow borehole network surrounding the site consisted of 12 4.5 Hz geophones with sampling rates of 500 Hz, located at various depths (0.2-1 km), enabled a better directivity analysis. We use this network to examine azimuthal differences observed in the spectral and temporal content of the earthquakes. We use the catalog by Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al. (2021), , who refined the original catalog by Kwiatek et al. (2019). Of the 5,456 relocated events, we analyze the 21 largest events ( 1.3-1.9) to investigate directivity effects.

Directivity Analysis
Directivity effects are identified by examining the source at different azimuths (e.g., Frankel et al., 1986). In the time domain, the observed source duration, that is, the width of the seismic moment rate (known as the apparent source time function, ASTF), is shorter at stations in the forward rupture direction. Alternatively, the corner frequency of the source spectrum, reflecting the observed high-frequency content, is inversely proportional to the source duration and is thus larger in the forward direction. Figure 2a illustrates the expected source duration, , and corner frequency, ci , at station for a unilateral earthquake for a given double-couple focal mechanism using the Haskell (1964) directivity model: 10.1029/2022JB025226 3 of 14 where is the earthquake rupture duration, is the rupture velocity, is the P-or S-wave velocity, and is the angle between the rupture and the ray path directions. The term cos can be expanded to cover the full focal sphere (e.g., Park & Ishii, 2015): cos = sin sin + cos cos cos( − ) where ( ) and ( , ) are the dip and azimuth angles of the rupture and take-off directions, respectively. Figure 2a shows that the forward rupture and slip direction, indicated by an arrow, experiences the smallest and largest ci . The then increases symmetrically along the focal sphere away from the rupture direction until it reaches its maximum at the backward direction, and vice versa for ci . Stations along the auxiliary plane observe constant and ci . In more complicated ruptures, such as bilateral ruptures, the resultant and ci patterns will be more complex and diverge from symmetry (e.g., Meng et al., 2021). However, any forward rupture directions will still experience smaller and larger ci than the rest of the focal sphere.
To analyze directivity effects through or ci , we first isolate the earthquake source component through deconvolution using a smaller, co-located Empirical Green's Function (EGF) earthquake (e.g., Bakun & Bufe, 1975;Holmgren et al., 2019;Mori & Frankell, 1990). Choosing the 21 largest events ( 1.3-1.9) as target earthquakes and focusing on the S-wave, we initially select EGF events for each target by examining the 50 closest earthquakes. Using a 0.6 s time window, we cross-correlate the waveforms of the EGF events with the target to ensure similar focal mechanisms and locations. We apply a 15-40 Hz, two-pole, two-pass Butterworth filter to each record, keeping any EGF seismogram with a cross-correlation coefficient (CC) ≥ 0.7 for future steps. Next, we obtain Fourier spectra using multitapering (Prieto et al., 2009) and compute the signal-to-noise ratio (SNR) for each record using a 0.6 s pre-P-phase segment as the noise. We consider frequencies between 7 and 212.5 Hz (85% of the Nyquist frequency). To select target-EGF event pairs, we require SNR ≥ 3 starting at maximum 10.5 Hz (1.5 × minimum considered frequency) and encompassing a minimum bandwidth of 50 Hz. With potential EGF events selected for each target event at each station, we begin our directivity analysis by examining station variations in ci .
The source spectrum of an earthquake can be described by the omega-squared model (Boatwright, 1980;Brune, 1970):   with magnitude proportional to the size of the circle. Red circles indicate the 21 largest events ( 1.3-1.9). The stimulation intervals are shown, color-coded based on stimulation stage (Kwiatek et al., 2019), and maximum horizontal stress ( max ) direction N110°E is indicated (Kakkuri & Chen, 1992 where Ω0 is the low-frequency plateau, is the frequency, is corner frequency, and is a constant defining the sharpness of the corner. We use the Boatwright (1980) model with = 2 . In the EGF method, the spectral ratio between the target and EGF events is: where subscripts and EGF indicate target and EGF. For each target-EGF record pair we compute the spectral ratio in the frequency domain. We discard any spectral ratios where the low-frequency ratio (Ω 0, ∕Ω0,EGF ) is less than 5.6, equivalent to a difference (Δ ) less than 0.5, to ensure the EGF is small enough. For each station with a minimum of four spectral ratios, we stack the normalized spectral ratios and use nonlinear least squares to fit the stack to Equation 4 and solve for (i.e., ci ). Following Viegas et al. (2010) and , the uncertainty bounds are defined as the minimum and maximum frequencies over which the variance is within 5% of the minimum variance. We also test a stricter Δ criterion of 1.0 (equivalent to a low-frequency ratio of 31.6), which is commonly used to ensure the delta function assumption holds (e.g., Huang et al., 2017;Onwuemeka et al., 2018). While this removed 5,471 individual spectral ratios out of 10,240 spectral ratios in total and several station stacks, we find that the final station results do not change significantly and stay within the uncertainty bounds from the 0.5 Δ results, thus deciding to keep the 0.5 Δ criterion. Of the 21 analyzed target earthquakes, 12 exhibit clear spatial variations in ci , spanning between 1.4 and 1.9. Three of these are shown in Figures 2b, 2d and 2f. See Figures S1-S21 in Supporting Information S1 for the analyzed events' waveform timeseries and spectral ratio stacks at each station.
We next investigate directivity effects in the time domain by deconvolving the targets by their EGFs to obtain ASTFs (e.g., Hartzell, 1978). First, to ensure the EGF duration is small enough to satisfy the delta function assumption and avoid underestimating the ASTF pulse width (e.g., Lanza et al., 1999), we use the and EGF outputs (Equation 4). Considering is inversely proportional to , we generate far-field STF pulses representing targets and EGFs of different magnitude and stress drop ranges and investigate how much smaller the EGF pulse width ( EGF ) is required to be to retrieve the target pulse width ( ) from the deconvolved ASTF. We find that ∕ EGF ≥ 2.5 resolves the ASTF pulse width to within 75% of the true target duration, thus requiring EGF∕ ≥ 2.5 to compute ASTFs. We also tried using a stricter criterion of EGF∕ ≥ 4.0 (resolving to within 85% of ) but found that the output ASTFs were close to identical albeit significantly fewer. Next, to ensure the ASTFs are stable, we use a water level correction of 0.001 (e.g., Mueller, 1985). To remove the high-frequency noise, we apply a 212.5 Hz low-pass filter (85% of the Nyquist frequency). Finally, for each target earthquake, we stack all normalized ASTFs for each station. We calculate values by estimating the width of the stacked ASTFs at 0.1 of the maximum amplitude (Courboulex et al., 2016). We define the uncertainty as the sampling rate (± 0.002 s).
Of the 21 target earthquakes, 10 events display directivity effects. These 10 are included in the 12 events with directivity from the ci results. The remaining earthquake with clear ci variation has too few resolvable station ASTFs to determine directivity. Station ASTFs are shown for three directivity events in Figures 2c, 2e and 2g. Based on visual inspection, eight of the 21 events also displayed complex ruptures: the ASTF shape deviates from a simple pulse (e.g., Figures 2c and 2g). Complexity in the frequency domain is difficult to model and typically appears as bumps altering the spectral shape (e.g., Holmgren et al., 2019;Uchide & Imanishi, 2016), resulting in biased ci estimates. However, while complexity is easier to identify in the time domain, we find fewer ASTFs in the Helsinki data set and thus the ci results provide overall better focal sphere coverage. Furthermore, identical trends for ci and suggests that ci is a useful indicator of forward rupture directions even if the rupture is complex. Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al. (2021) determined focal mechanisms for 191 earthquakes from station polarities using a cross-correlation technique and manual inspection. They identified three distinct families with oblique reverse faulting as the dominant source mechanism type, a mechanism not easily resolved with the available station coverage. Due to the combination of deep microseismicity (∼6.1 km) and relatively small epicentral distances to stations (≤ 9 km), the focal sphere station coverage is clustered toward the center of a compressional quadrant. Thus, constraining the locations and orientations of the nodal planes depends on the few stations at the network edges with lower SNR. Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al. (2021) found that Family 1 was fairly well constrained with consistent polarities, but the remaining two families displayed more polarity variation between events and were less stable.

Focal Mechanism Reassessment
We refine the focal mechanisms using a Principle Component Analysis (PCA) approach (Vavrycuk et al., 2017), which combines classical amplitude and waveform-based moment tensor (MT) inversion techniques and is suitable for microseismic, high-frequency events with low SNR. First, we apply a 5-20 Hz, two-pole Butterworth filter to each target earthquake displacement record at each station, extracting 2.0 s time windows centered on the P-wave arrival. For each event, we taper the station waveforms, align them using cross-correlation, and extract the common source wavelet (see details in Vavrycuk et al., 2017). Finally, we retrieve MTs through classical MT-amplitude inversion (hybridMT, Kwiatek et al., 2016), using the PCA coefficients as input. We calculate stable double-couple constrained MTs with small non-double couple components for 19 of the 21 events (see Figure 3 and Table S1 in Supporting Information S1). The additional usage of wavelet-based amplitude input leads to more homogeneous oblique reverse focal mechanisms consistent with Family 1 of Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al. (2021) (see Figure S22 in Supporting Information S1). The changes in focal mechanisms are primarily due to improvement of ambiguous P-wave polarity from stations located close to nodal planes. The obtained mechanisms, regardless of nodal plane, are well oriented within the stress field (cf. Figure 10 in Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al., 2021). We find that the PCA approach is well-suited for MT inversion of small earthquakes by building on amplitude-based inversion and yet being less complex than a full-waveform approach (see discussion in Bentz et al., 2018). For example, because it relies on the area under the P-wave pulse, it is less sensitive to directivity patterns than the classical amplitude-based inversions. Furthermore, path and site effects are diminished because the sensors are located in boreholes drilled into crystalline Precambrian Svecofennian basement rock with very low attenuation (e.g., Kwiatek et al., 2019), allowing for reliable focal mechanism solutions even though the earthquakes' small magnitudes require source analysis in the higher frequencies (here, we focus on the frequency range 5-20 Hz). Eulenfeld et al. (2022) analyzed attenuation effects using the Helsinki induced events and found that the quality factor is larger than 1000 for frequencies above 10 Hz. Abercrombie (1995) found that using deep borehole sensors removed the severe attenuation that occurs in the upper kilometers in a granite batholith at the Cajon Pass, California. We note that Ide et al. (2003) reported strong path and site effects in a deep borehole at the Long Valley caldera in California, however, the attenuation properties of a volcanic caldera are significantly different to a crystalline batholith.

Combining Directivity and Focal Mechanisms
Next, we combine the focal mechanisms and observed directivity to identify fault planes from auxiliary planes and invert for the earthquake rupture directions using Equation 1 (e.g., Jost et al., 1998;Li et al., 1995). Because the observed directivity has better station coverage in the frequency domain than in the time domain, we use the observed ci estimates converted to instead of the observed estimates directly. A conversion commonly assumed in the literature is simply ≃ 1∕ (e.g., Archuleta & Ji, 2016;Van Houtte and Denolle, 2018). Empirically, Prieto et al. (2017) found a similar conversion ( = 0.94∕ ) for a 75 km deep 4.8 earthquake. Meanwhile, Hisada (2000) theoretically determined = 1∕(2 ) , while Denolle and Shearer (2016) Tomic et al. (2009) computed both spectral ratios and ASTFs using the EGF method to investigate stress drops of ≤ 2.1 induced earthquakes in Brazil. After converting their reported averaged rupture radii back to and finding the mean duration for each event (see Tables 2 and 3, Tomic et al., 2009), their earthquakes also show a = 1∕( ) relation, where varies between 1.06 and 1.58 between the events. Here, we empirically investigate the relationship between and to find the most suitable conversion. Using all the event-station pair stacks with both measured ci and estimates, we determine a linear empirical relationship between the values (see Figure 4a). To ensure the values are not too small to be resolvable considering the 500 Hz sampling rate, we first synthetically examine the minimum number of samples required inside a pulse to estimate the pulse width (see Figures S34-S38 in Supporting Information S1). For a symmetrical pulse, we find ≥ six samples can resolve the pulse width within 5% of the true width. For the Helsinki data, which has a sampling rate of 500 Hz, a minimum of six samples corresponds to a minimum ASTF pulse width = 0.012 s. This minimum limit is indicated as a dashed horizontal line in Figure 4a and we find that the versus 1∕ ci relationship below this limit does not visually follow the same trend as the remaining data. We find the following relationship between the measured and ci estimates:  Denolle and Shearer (2016). To visually inspect this relationship further, we also include the event-station ASTF and spectral ratio stacks in Figures 4b and 4c with their and ci estimates shown (dark diamonds). Using the corresponding ci or measurements, we calculate the resultant or ci using Equation 5 (pink squares). We also include the resultant or ci estimates using the commonly assumed ≃ 1∕ relationship (purple circles). As can be seen, the ≃ 1∕ conversion leads to overestimation in both and ci at longer durations and lower corner frequencies, respectively.
For the inversion, we constrain the possible rupture directions to the two nodal planes, assuming unilateral rupture. In short, we find the optimum rupture direction for each nodal plane and then determine the earthquake's optimum rupture direction by comparing the root-mean-square (rms) residuals of the two directions. Figure 5 illustrates the inversion process for the largest event, 6389 ( 1.9). We let the strikes and dips obtained from the focal mechanism solutions constrain the possible rupture directions in the inversion and solve for a best-fit rupture angle on the planes for the rupture direction to retrieve from Equation 1. This is done by combining the rupture angle, strike, and dip to give us a plunge and trend of the rupture direction on the plane. The plunge and trend are then used to calculate the distance along the sphere to each station, that is, . For each of the two nodal planes, we solve for the rupture angle, , and using the nonlinear least-squares inversion in MATLAB, also extracting the Jacobian to retrieve the 68% confidence interval for each parameter as an uncertainty measure. Thus, our objective function becomes: where values are obtained by converting the observed ci using Equation 5. After the two nodal planes' optimum rupture directions and their uncertainties have been found, we conduct a grid search around their rupture-angle uncertainty bands to evaluate which rupture direction is the earthquake's most-likely rupture direction (and thus which nodal plane is the fault plane). For each evaluated rupture angle in the grid search (see Figure 5e), we let and be free parameters and redo the inversion to compute the rms of the residuals. If one of the planes has lower rms(residuals) over its full confidence interval, we select this as the event fault plane and rupture direction. If the rms(residuals) overlap over the confidence intervals, both nodal planes and rupture directions are kept (see Figures S23-S33 in Supporting Information S1 for more event examples).
Out of the 12 events with directivity, we constrain one best-fit rupture direction for 10 events and two possible rupture directions for two events (see Table S2 in Supporting Information S1). The inverted rupture direction(s) for each earthquake is shown in Figures 6a and 6b, where the rupture direction is indicated and colored based on angle with respect to the well. Most of the events appear to rupture parallel to the well. Two events out of the 10 events with one resolvable rupture directions rupture away from the well: events 5695 ( 1.6) and 6389 ( 1.9, the largest event) both occurring during stimulation stage 4. However, when considering the ambiguity in nodal planes and the 52-m catalog location uncertainty (Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al., 2021), the rupture direction with respect to the well is more uncertain for event 5695, which occurred close to the well (44 m away). To investigate this further for each event, we let the rupture direction stay the same but let the earthquake location change along a sphere centered at the event hypocenter with a 52-m radius. We compute the rupture angle at each point on the sphere to examine how much it changes considering the location uncertainty (see Figure S39 in Supporting Information S1). Unsurprisingly, we find that the events within 55 m of the well (events 4364, 5695, 5386, and 7641) are much more sensitive to the location uncertainty and their rupture angles are less stable. However, the largest event, event 6389, occurred 132 m away and thus still unequivocally ruptured away from the well, not toward the well as found at other geothermal sites (e.g., Folesky et al., 2016;Kiraly-Proag et al., 2019). We also examine the evolution of the rupture angle with respect to time and distance from the well (see Figure 6c), finding no clear trend.
Finally, we use the rupture angles on the fault planes to investigate the prevalence of mode II (in-plane shearing, slip and rupture are parallel) and mode III (out-of-plane shearing, slip and rupture are perpendicular) failures among the earthquakes. For each event, we find the angle between the rupture angle and the rake angle from the focal mechanism solutions (Figure 7a), where 0° and 180° (i.e., rupture parallel to the slip direction) indicates ). (c, d) Same as (a, b) but for nodal plane 2 for the same event. (e) Grid search around the optimum rupture angle direction's 68% confidence interval plotted against the root-mean-square (rms) residual. For this event, plane 2 has the lowest rms(residual) and is the preferred nodal plane.  Figure 7b shows a polar histogram of the angles found, indicating a preference toward mode II failures among the 12 directivity earthquakes.

Discussion
Knowledge about the fault plane orientations and rupture directions of fluid-induced seismicity offer important insight into the driving forces governing earthquake behavior. Such information is essential to estimating key parameters such as maximum rupture size and : rupture growth models have shown that pore-pressure driven seismicity is more likely to be constrained within the perturbed zone (e.g., Dempsey & Suckale, 2016;Galis et al., 2017;Gischig, 2015;McClure & Horne, 2011;Shapiro et al., 2011), while low-pressure and high fault-stress locations are more probable to host runaway ruptures (e.g., Galis et al., 2017;Gischig, 2015). One distinct earthquake behavior related to the governing force is the rupture direction with respect to the injection wells, where rupture toward the well is more probable in pore-pressure driven environments with high injection pressures (Dempsey & Suckale, 2016). In observational studies, this has been observed particularly in the largest events at both deep geothermal and high-rate wastewater injection sites ( > 1.8, Folesky et al., 2016;5.1-5.7, Lui & Huang, 2019; 3.3, Kiraly-Proad et al., 2019). The induced seismicity from the Helsinki geothermal stimulation is primarily pore-pressure driven, as evidenced by max being relatively low and following the Galis et al. (2017) model for stable, pressure-controlled injection and the seismicity mainly constrained to the diffusion cloud (Kwiatek et al., 2019). However, we do not observe a directivity bias for rupture toward the well amongst the largest Helsinki events. Indeed, only one event ruptures toward the well (event 4364, 1.7, see Figure 6), although when the proximity to the well (55 m) and the catalog location uncertainty (52 m; Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al., 2021) are considered, the rupture angle becomes less stable (see Figure S39 in Supporting Information S1). If we also consider the two events with two possible rupture directions, one of the events (6246, 1.5) has a plane with a preferred rupture toward the well that is more stable. Thus, at best, this results in two out of 12 events which appear to rupture back toward the well. Dempsey and Suckale (2016) examined directivity bias as a distribution of directivity observations from different pressure evolution scenarios, also determining the smallest number of samples required to establish directivity bias at a site using Kolmogorov-Smirnov testing. They found that for the scenario applicable to EGS sites

(b)
("advancing-front" scenario where permeability enhancement dominates), the directivity bias was so strong that only five observations were required to detect it (see Figure 14 in Dempsey & Suckale, 2016). Thus, the Helsinki site appears to not exhibit the expected rupture direction behavior, especially considering that the same number of events rupture away from the well (events 6389 and 5695) as toward. Instead, the preferred rupture direction seems to be roughly parallel to the well, displayed by eight out of the 12 events.
One common assumption in rupture growth models is that earthquake nucleation initiates at the pore pressure front (e.g., Dempsey & Suckale, 2016;McClure & Horne, 2011), which has a steeper gradient at the beginning of stimulation. Meanwhile, the studied Helsinki events occur at the earliest during stage 2, most occur during stages 4-5, and they locate toward the bottom of the seismicity cluster ( Figure 6). Thus, this suggests these events occurred in a highly stressed rock volume continuously supplied with new fluids. In the absence of a distinct pore pressure front, it is possible that the pressure gradient no longer controlled the directivity of the largest events to the same extent as during early stimulation stages. Instead, later events could be expressing the local stress heterogeneities, stress concentrations, and structure of the fracture zone, allowing for alternative rupture directions. Unfortunately, only one event (event 1585) from stimulation stage 1 was large enough to be included in this study and it did not display any clear directivity, preventing us from examining any possible rupture directivity bias during the early stages when the pore pressure gradient was steeper. For the later stages, the Helsinki earthquakes show a preference toward parallel directivity with respect to the well, possibly indicating that the local fracture zone is controlling the rupture process. Interestingly, this is the opposite behavior to what Folesky et al. (2016) observed with their smaller events. While their largest events ( ≥ 1.8 ) ruptured back toward the well, the smaller events displayed an increasingly clearer trend of orientation with distance away from the open hole. They suggested that the increased pore pressure controls the rupture process in the vicinity of the borehole (<100 m), resulting in no clear trend in rupture directions, and that either the stress field or geological structure takes over once the pore pressure influence decreases. For the Helsinki events, seven out of the nine events within 100 m (and the increased pore pressure cloud) display parallel directivity to the well (see Figure 6c).
Rupture growth models also typically consider a single, isolated fault or fracture that allows for clear stress concentration and transfer (e.g., Dempsey & Suckale, 2016;McClure & Horne, 2011), rather than a fracture network in which fluid can take more complex paths, especially toward the end of stimulation. Indeed, in Helsinki, the refined focal mechanisms obtained using the PCA approach indicate reactivation of a network of parallel fractures, supporting earlier interpretations (Hillers et al., 2020;Kwiatek et al., 2019;Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al., 2021). Hydraulically stimulated fractures are generally a combination of hydrofracking (mode I, tensile opening) and hydroshearing (mode II and III, in-plane and out-of-plane shearing, respectively) (e.g., Frash et al., 2019). Reactivation of a pre-existing fracture network, which is a likely mechanism of the seismicity during the 2018 Helsinki simulation (cf. discussion in Kwiatek et al., 2019;Leonhardt, Kwiatek, Martínez-Garzón, Bohnhoff, et al., 2021), tends to predominantly result in hydroshearing as it allows for slip at lower pressures than hydrofracking (e.g., Pogacnik et al., 2016). This is supported by the small non-double couple components observed in our moment tensor solutions. Moreover, we find that the Helsinki microseismic events are predominantly mode II failures rather than mode III (Figure 7). Generally, strike-slip faults are initially modeled as mode III and thrust/normal faults as mode II (e.g., Gudmundsson, 2014;Oglesby et al., 2000). Thus, a preference toward mode II failures agrees with the oblique reverse faulting source mechanisms observed at the Helsinki site. In contrast, Dempsey and Suckale (2016) modeled strike-slip faults in their numerical analysis, assuming mode III failures. We also investigate if we would be able to detect a mode III rupture given the dominant focal mechanism and station configuration in Helsinki by modeling the expected directivity patterns from both pure mode II and mode III failures and retrieving the range in found for each mode (see Figure S40 in Supporting Information S1). We find that for mode II ruptures we will have a coverage of 78°, while mode III produces a coverage of 55°, thus it easier to detect mode II failures. However, the event closest to a mode III failure (event 7641, see Figure 7) also had the smallest change in observed ci over its stations (22-37 Hz) with a low resolved rupture velocity ( = 0.20 , see Table S2 in Supporting Information S1) and was still resolvable in terms of directivity and rupture direction (see Figure S31 in Supporting Information S1). Radiation pattern is another factor that could affect the resolvability of mode III failures. Any stations located in the direction of the B-axis (or null-axis) would observe damped amplitudes due to unfavorable radiation pattern of the S-wave and could potentially bias the directivity pattern. However, because the Helsinki events do not have any station coverage over the B-axis, we do not consider this to significantly affect our results. Investigating the implications of a mode II preference over mode III is outside of the scope of this article, instead these results can be used by geomechanical modelers as observations to better understand fracture mechanics of fluid injected into fracture zones.
Finally, in the directivity inversion, we generally obtain low (∼ 0.19-0.66 ) rupture velocities (see Table S2 in Supporting Information S1). Low may indicate that more elastic strain energy is spent on fracture reactivation rather than radiation of seismic waves (e.g., Kanamori & Rivera, 2004;Kwiatek et al., 2011;Prieto et al., 2017), which seems feasible in a complex distributed fracture network. However, can be difficult to constrain in directivity analyses. For example, Abercrombie, Poli, and Bannister (2017) found that resolution and frequency limitations significantly affect . While we synthetically investigated the minimum resolution for a symmetrical pulse (see Figures S34-S38 in Supporting Information S1), the ASTF pulses could vary in shape and might in some cases be under-or overestimated. is also sensitive to the assumed ∝ 1∕ relationship. If we use = 1∕ instead of Equation 5, the rupture directions of the nodal planes stay the same but the obtained range increases to 0.37-0.90 , which results in the median increasing from 0.40 using Equation 5 to 0.63 . While we deemed = 1∕ not appropriate for the Helsinki earthquakes due to overestimation at longer and lower (see Figure 4), our ASTFs pulses were measured at 0.1 times the maximum amplitude due to the noise levels on either side of the pulse. For a symmetrical pulse, this results in a 5% smaller width (see ). However, considering the uncertainties were set to ± 0.002 s (i.e., the sampling step) and a 5% mismatch for the widest observed pulse (0.033 s) would lead to 0.0016 s error, it falls within the uncertainty. Finally, non-unilateral ruptures may lead to underestimated rupture velocities (Tomic et al., 2009). While we observe complex station ASTFs for eight of the events, indicating an additional rupture of a possible subevent or asperity that could lead to biased ci estimates, we did not identify clear bilateral ruptures.

Conclusion
We investigated the rupture behavior of fluid-induced microseismicity from the 2018 St1 Deep Heat geothermal hydraulic stimulation campaign in Helsinki using directivity and focal mechanisms, resolving rupture directions for nine of the largest events in a pore-pressure driven environment. We find rupture directions are variable with respect to the injection well, with one event rupturing toward, two rupturing away (including the largest event 1.9), and the remaining rupturing roughly parallel to the well. This is contrary to predictions by rupture growth models, which commonly assume a strong pressure gradient exists. These events occurred toward the end of stimulation and near the well, indicating that rupture directions may be more random when the pore-pressure gradient is weaker due to prolonged injection. Finally, we find that the events exhibit a preference toward mode II failure (parallel slip and rupture directions). These observations provide new constraints for mechanical models of rupture growth in pore-pressure dominated settings.

Data Availability Statement
All waveform data used is available in the associated data publication (Holmgren et al., 2022). Full earthquake catalog can be obtained through .