Induced seismicity due to hydraulic fracturing near Blackpool, UK: source modeling and event detection

Monitoring small magnitude induced seismicity requires a dense network of seismic stations and high-quality recordings in order to precisely determine events’ hypocentral parameters and mechanisms. However, microseismicity (e.g. swarm activity) can also occur in an area where a dense network is unavailable and recordings are limited to a few seismic stations at the surface. In this case, using advanced event detection techniques such as template matching can help to detect small magnitude shallow seismic events and give insights about the ongoing process at the subsurface giving rise to microseismicity. In this paper, we study shallow microseismic events caused by hydrofracking of the PNR-2 well near Blackpool, UK, in 2019 using recordings of a seismic network which was not designed to detect and locate such small events. By utilizing a sparse network of surface stations, small seismic events are detected using template matching technique. In addition, we apply a full-waveform moment tensor inversion to study the focal mechanisms of larger events (ML > 1) and used the double-difference location technique for events with high-quality and similar waveforms to obtain accurate relative locations. During the stimulation period, temporal changes in event detection rate were in agreement with injection times. Focal mechanisms of the events with high-quality recordings at multiple stations indicate a strike-slip mechanism, while a cross-section of 34 relocated events matches the dip angle of the active fault.

through temporary or long-term suspension of production, especially in regions that are not seismically active. For instance, prior to the seismicity studied herein, magnitude 1.5 and 2.3 earthquakes occurred near Blackpool, UK, during a previous shale gas reservoir development in 2011 (British Geological Survey 2017). These were the largest events in a sequence of induced tremors and led to a temporary UK-wide shutdown of hydraulic fracturing activities until 2018. While the majority of induced seismicity is minor, in rare cases events can be sufficiently large to pose a risk to structures at the surface (Grigoli et al. 2018). For example, a M W 4.6 earthquake was caused by fluid injection during hydraulic fracturing in northeast British Columbia in 2015 (Clay 2015).
Monitoring induced seismicity is, therefore, clearly very important and crucial to the safe and economical operation of any site where induced events are expected, whether shale gas, geothermal or carbon capture storage López-Comino et al. 2017;Karamzadeh et al. 2019). The distribution of event locations, and their mechanisms can lead to the identification of hidden faults intersecting or close to the injection well-bores. In addition, the source characteristics of the events provides a criterion to discriminate fault related and fracture related events (Kratz et al. 2012). The occurrence of larger microseismic events is an indication of fault reactivation (Davies et al. 2013) due to the interaction between injected fluid and the fault, either by permeation of the fluid into the fault plane, or perturbation of the stress state of the fault. Furthermore, microseismic monitoring allows us to infer the fracture propagation length and the extent of stress changes (Wilson et al. 2018).
Generally, recording and locating weak seismic events require high-quality seismic data acquisition with a very dense network of seismic sensors at surface and/or depth, which is very often not available or restricted to the operator. Besides data availability, the hydrofracking-related activity (such as traffic) increases seismic background noise level and environmental vibrations Lopez-Comino et al. 2018). This makes seismic monitoring from the surface more challenging, as small magnitude earthquakes might be buried in the background noise and not be detectable using techniques based on time-domain amplitude inspection. Template matching (TM) techniques offer a solution to detect such weak seismic events (Gibbons and Ringdal 2006;Meng et al. 2018) and have been shown to increase detection capabilities in other applications (Yoon et al. 2015;Beaucé et al. 2019;Skoumal et al. 2016).
In this study, we aim to show the capabilities of surface stations alone to detect and characterize the induced seismicity in Blackpool during hydraulic fracturing of the second well PNR-2, using modern waveform-based analysis tools. By utilizing a sparse network of surface stations, a number of seismic events are detected using template matching technique. The temporal behavior of the detections are studied and compared with the stages of hydraulic fracking operations taking place at about 2-km depth. In addition, we applied a full-waveform moment tensor inversion to study the focal mechanisms of larger events (M L >1) and used the double-difference location technique for events with high-quality and similar waveforms to obtain accurate relative locations.

Induced seismicity near Blackpool
As a consequence of hydraulic fracturing activities in the subsurface near Blackpool, several sequences of earthquakes occurred and were reported by the BGS. In 2011, on April 1 and May 27, two earthquakes, with magnitudes 2.3 and 1.5 M L , respectively, were detected and felt in the Blackpool area. These events were part of a sequence that occurred between March 31 and May 27 comprising a total of 52 tremors, with magnitudes ranging between M L −2 and M L 2.3 (Clarke et al. 2014). On October 18, 2018, fracking restarted at the nearby Preston New Road site, after a 7-year UK-wide moratorium, and consequently earthquakes with magnitude of M L 1.1 and M L 1.5 occurred on October 29 and December 11, 2018. These two events were the largest and felt earthquakes in a series of 57 events with magnitude between M L −0.9 and M L 1.5 (Galloway 2019). After restarting operations in 2019, at a second well at the same site, the occurrence of two earthquakes with M L 1.6 and 2.1, in a sequence of 137 events, led hydraulic fracturing at the site to be suspended on August 24. The largest event, M L 2.9, occurred 2 days later on August 26. The largest event was associated with reports of minor damage to local belongings and was assigned EMS-98 intensity VI by the BGS. (Edwards et al. 2021) found that this was likely an overestimate, with intensity V being more representative of the reported and modelled effects. As a result of a regulator-led scientific review of induced seismicity near Blackpool, the UK government again imposed an indefinite UK-wide suspension to onshore hydraulic fracturing in late 2019.

Hypothesis about induced seismicity near Blackpool
Analysis of the sequence of seismic tremors in 2011, using 4 local stations at 1-to about 2-km distance from the epicenters, revealed the location of a hidden fault that was activated during hydraulic fracturing (Clarke et al. 2014). Source modeling of the biggest event in the sequence indicated a strike-slip mechanism with left-lateral motion (Clarke et al. 2014).
Induced seismicity related to the 2018 sequence was recorded by a dense surface network of 26 surface broadband and short-period sensors and 24 downhole geophones (15 Hz) at the depth at which injection took place. Seismicity was strongly clustered in space and time, and associated with known periods of injection. In addition, the observed seismicity migrated from west to east, in agreement with the spatial locations of different stages of hydraulic fracturing (Baptie 2018). From analysis of six of the largest events of 2018 stimulation, Clarke et al. (2019) obtained focal mechanisms implying either left-lateral strike-slip on a near-vertical fault striking northeast-southwest, or right-lateral strike slip on a near-vertical fault striking northwest-southeast. Based on stress analysis, they concluded that the maximum and minimum horizontal stresses are oriented in north-south and east-west directions, respectively; thus, the main plane causing the events is optimally oriented for left-lateral strikeslip motion. Analysis of induced seismicity related to the 2019 hydrofracking stimulation by (Kettlety et al. 2020) concluded that the aftershocks of the M L 2.9 event were distributed along a near vertical plane with a strike and dip of 130 • and 80 • , respectively. In addition, the focal mechanism of the mainshock determined from surface array data using first motion polarities of both P-and S-waves and P-to-S-wave amplitude ratios indicated a strike of 127 • , dip of 84 • and rake of −160 • , i.e, the orientation of the northwest-southeast nodal plane consistent with the plane fitted to the aftershocks. According to Kettlety et al. (2020), the spatial distribution of events in 2018 and 2019 indicates that the seismicity occurred on two different faults, as the stimulated well in 2019 was 100 m above the previous one. The fault reactivated in 2018 had a northeast strike, whereas the fault reactivated in 2019 had a southeast strike (Kettlety et al. 2020).

Processing method
In this section, first the analysed data related to the 2019 microseismic sequence are introduced, and then the employed processing methods including event detection using template matching technique, source modeling technique and event location are discussed in detail. For each processing method, the result of application on the data is also given.

Data and data quality
The data includes continuous seismograms recorded at stations of two local broadband seismic networks, 'UK Array', with network code UR (Bristol University 2015), operated by the British Geological Survey (BGS) and 'Seismic Data in Northwest UK', operated by the University of Liverpool, network code LV (Liverpool University 2014). Figure 1 shows the location of the nine selected stations used for event detection and focal mechanism modeling. The sampling rate of most stations in the UR network is 200 sps; however, stations AQ07 and AQ09 are exceptions with 100 sps. The sampling rate of all LV stations is 100 sps. The data, spanning about 7-week duration (from July 30, 2019, to Sep 20, 2019), is processed, starting from 3 weeks before the first 'large' event in the sequence (M L 1.6 on August 21) until 3 weeks after the biggest event, M L 2.9 (August 26). While data of all stations are explored to detect events, we mainly focus on the nearest station to the observed seismicity, site AQ04. In the processing time period, BGS reported 133 events with magnitude above zero within 1-km radius around the location of the M L 2.9 event. Locations of the earthquakes are shown by black stars in Fig. 1. The red inlay in this figure includes the epicenter of the events relative to the largest event (T5) which is shown with a red circle. The number of phases which contributed to the BGS location procedure is between 4 and 28 phases from recordings of UR and the operator networks. Time, magnitude and location of the 6 largest events in the time period investigated in this study are listed in Table 1. The focal mechanisms shown are calculated for each of template events except T6 which lacks of enough high-quality recordings (see Table 1) The probabilistic power spectral densities (PPSD) (McNamara and Buland 2004) of seismic recordings indicate the quality of signals. For all stations, modal values of PPSD are shown in Fig. 2. We note that the stations of LV network were unusable between August 8 and August 17 due to segmented data and several gaps. In general, the PPSD mode curves are rather close to the high noise level for the processed time duration of all stations. The difference in noise level between stations can reach up to 30 dB, while AQ04 experiences the lowest noise level at high frequency band (above 5 Hz). Stations of the LV network show more fluctuations and exhibit the highest level of noise in high frequency band. The high and variable noise levels are due to the predominantly urban or residential and near-coastal localities of the seismic networks, which for cases where induced seismicity poses a risk, will be unavoidable.

Event detection using template matching
Hydrofracking sites are often noisy and the induced events are of relatively low magnitude. Most of the events are, therefore, not detectable at the surface using a detection technique based on single-station time-domain amplitude inspection. However, template matching techniques offer a possible solution to overcome this high noise situation, increasing the detection capability of weak events at the surface (Gibbons and Ringdal 2006;Meng et al. 2018). Template matching methods make use of the waveforms of already extracted events as templates, e.g. a few larger events with a sufficient quantity of high-quality phase picks. These templates are used to scan through the continuous waveform data, prior or posterior to the template timestamp, and to find new similar events. The comparison is made by means of correlation measurement, and detection is declared once the value of the correlation between a template waveform and the piece of signal under analysis exceeds a predefined threshold. The processing steps of the algorithm are briefly introduced in Table 2. This technique is particularly useful when the noise level is high compared to the size of the events of interest, and the seismicity happens in a relatively small seismogenic zone, which is the case for induced seismicity related to fluid injection (Gibbons and Ringdal 2006;Shelly et al. 2007).
By applying TM approach to our data set, we anticipate that weak events which might be buried in seismic background noise are detectable and we can extract as many events as possible to compile an event catalog with lower detection threshold. The aim is to see if we can detect and process the smaller events identified by few stations at surface and link observations at surface with the actual process that is happening at seismogenic depth. We assume that the events are generated in a region smaller compared to the sourcereceiver distances so the time difference between Sand P-phase arrivals are consistent for all events at one station. We used EQcorrscan (Chamberlain et al. 2018) which is an open-source software package written in Python, for detection and analysis of seismicity, using an eight-core processor computer. The main processing steps are as below: -The TM algorithm is applied to recordings of the nearest station to the observed seismicity, AQ04. -All detected events are clustered using a waveform-based clustering approach, and a number of trusted events are extracted. -TM algorithm is applied to other stations and the detectability of the small events from distant stations is studied. -For each detected event, the P-phase onset time is re-picked. -Events with multi-station picks are extracted for further analysis.
It should be noted that while event detection is performed by measuring the stack of cross-correlation of the 3-component traces, P-phase re-picking is performed simply by calculating the correlation of each component of a detected event with the same component of the template event, and measuring the time shift. As a result, for a detected event, there might be no P-phase picking in one component if the value of the correlation on that component is not above the threshold value of 0.4 (out of 1). If all 3 components of a detected event correlate well (greater than 0.4) with the related template waveforms, then 3 identical P-phase picking times are measured.

Templates
The six larger events, with magnitude greater than M L 1, are used as 3-component templates at each station. The time window of a template signal starts 0.1 s before P-phase onset time and includes both P-and Sphases signals, e.g. the window length for the station  Step 1: Making templates Inputs: 1: Continuous waveform (3-C) of a single station 2: Catalog of the template events (M L > 1) including P-and S-phases arrival times, t i p and t i s , i: template index 3: Processing parameters filter band, template length (t L ; t i L > t i s − t i p ) Process: 1: Pre-processing of the input waveforms filtering, mean removal 2: Cutting template waveforms from t i Template waveforms T i, i = 1 : n, n: number of templates Step 2: Template matching Inputs: 1: Continuous waveforms (CW) 2: Template waveforms 3: Processing parameters same as parameters used in Step 1 4: Threshold type and value (T hr) Process: 1: Pre-processing of the input waveforms 2: Computing cross-correlation of templates and CW CC(t) i , t: time instance of CW, i: template index 3: If CC(t) i > T hr then declare a detection 4: Keeping only the detection with higher CC in case of simultaneous detections with multiple templates Output: List of detected events and corresponding templates AQ04 is 2.0 s. Waveforms of template events at station AQ04 are plotted in Fig. 3, applying a band-pass filter of 20-80 Hz. Waveforms are aligned in time, based on the arrival time of P-phases and are scaled to the maximum value of each trace. According to Fig. 3, scattered coda phases coming after the P-phase of T5 (the M L 2.9 event, Table 1) show higher amplitude than the P-phase on the vertical component, and this is obviously different from the other templates. Similarly the S-phase coda shows higher amplitude than the first arrival S-phase on the horizontal components.
In addition, strong P-phase coda waves are visible on vertical component traces of T1 and T3, whereas such phases are relatively weak on T2, T4 and T6 (Table 1). Spectral properties of the templates are studied over a frequency range of 1 to 100 Hz using spectrograms calculated by Continuous Wavelet Transform (CWT) (Stark 2005). This analysis further refines the appropriate frequency bands to investigate using TM technique. Figure 4 shows the 3-component traces of templates T1, T5 and T6 and corresponding CWT coefficients. The traces are normalized so that the maximum value of each trace is 1; however, CWTs are calculated on unnormalized traces (band-pass filtered between 1 and 100 Hz). Comparing spectral behaviour of the east-west (E) components, we see that T5 shows higher but narrower frequency content of S-phases (in about 1 s) compared to the other two events. The second distinctive feature is observed on the northsouth horizontal component (N) when compared to the other templates: template T5 shows much lower frequency energy in almost the entire event time span, with a maximum value at about the onset of S-phase. The predominant frequency of the S-phase (i.e. frequency components with higher wavelet coefficients) is slightly smaller for T6 compared to T1. In vertical components (Z), the highest signal energy of T1  High frequency (20-80 Hz) and T3 are related to the S-phases, whereas for T5 the P-phase carries the most energy.
Waveform similarity between all template pairs are investigated by measuring the cross-channel correlation between two frequency bands. We then interpret the outcomes of different frequency bands. The values of cross-channel correlations are listed in Table 3. In the low frequency band (2-20 Hz), templates T1, T2, T3, T4 are highly correlated to each other, i.e. cross-channel correlation values are greater than 0.8, while in the high frequency band (20-80 Hz) template T3 shows less similarity to the others. Templates T5 and T6 are not correlated to the other templates or to each other, in both frequency bands. Template T5 shows slightly higher correlation to the T1-T4 and T6 in the high frequency band (correlations coefficient 0.31-0.37) than the low frequency band (0.2-0.28). On the other hand, template T6 shows higher correlation to T1-T5 in the lower frequency band (0.25-0.56) compared to the higher frequency band (0.0-0.17).

Template matching (AQ04)
Median absolute deviation (MAD) of daily stacked 3component cross-correlation traces is used as a threshold parameter in the TM detection procedure. The TM method is applied employing different threshold values and band pass filters, some of the outcomes are listed in Table 4, where the total number of detections, detection thresholds, frequency bands and the number of common detections with BGS catalog are provided. This table highlights that the frequency band of a template event influences the number of detections, as the value of cross-correlation between events is frequency dependent. The frequency band of 20-80 Hz is an appropriate frequency band for event detection and the threshold value is selected to be 9 MAD in the detection procedure, since using these values we obtain a smaller number of missing events compared to BGS catalog. By using 20-80 Hz filter band, 2217 events are detected (see Table 4). Including lower 2-20; 20-50 9 2956; 1501 20 (15%); 9 (0.06%) frequency content to the event detection procedure leads to a decrease in number of detections. For instance using band pass filter of 2-20 Hz and 10-80 Hz, the total number of detections is 372 and 578 events respectively. In addition, in the low frequency band of 2-20 Hz T5 does not make any detection. In Fig. 5, the cumulative number of detections at station AQ04 for each template event is shown for a frequency band of 20-80 Hz. The gray area in this figure shows the time period when the fluid injection was activated, between August 15 and August 23. The sum of detections made by templates T1, T2, T3 and T4 are also plotted in Fig. 5 as according to Table 3 those templates have waveform similarities. Temporal jump-like changes to the number of detections are observed using all templates. The lowest and highest number of detections are related to T2 and T3, respectively. During hydraulic fracturing (inside the gray area), the growth rate of detections with individual templates as well as overall cumulative curve exhibit fast changes in some periods, which are in agreement with the timing of pumping stages (dashed vertical lines). It is interesting that the general pattern of detection rate for T5 and T6 is more similar to one another than to the other templates. Templates, T1, T2, T3 and T5 show similar patterns of detection rate before the occurrence of T1. The number of detections with T5 and T6 increases during the hydraulic fracturing until the occurrence of T1, after which, for both cases, the cumulative number of detections is almost flat. In addition, slight changes are observable on all curves at the occurrence time of template events. For example, after the occurrence of T1, the number of detections with T1 (green curve) increases.
Results of event detections at station AQ04 are summarized in Fig. 6. September 6 was excluded from the analysis because of a data gap. The highest number of detections in a single day (128 detections) happened on August 20, 1 day before the magnitude 1.6 event (T1), and at which point the sequence of events with magnitude above M L 1 starts. A temporal jump in cumulative detection also occurs at that day for each template event and overall detections (see Fig. 5). Figure 6b shows the value of 3-component cross-correlation, i.e. stacked trace, at a detection time against the time of each detection. The total number of detections in the analysed time period is 2217 events, 729 of which (30%) happened during the night hours. Figure 6b also indicates 130 events that are reported in the BGS catalog.  Results of applying template matching event detection technique on seismic recordings at stations AQ04. In (b) each circle shows a detected event and the gray and white stripes distinguish between night (19:00 to 7:00 UTC) and day (7:00 to 19:00 UTC) hours

Clustering of detected events
Waveforms of all 2217 detections are further explored to extract similar events. Detected events are clustered by setting the average cross-correlation (cc) threshold to 0.6. The total number of re-clustered events is 485 extracted from all detections. By checking waveforms of all extracted clustered events manually, a cluster of 43 detections was recognised to be identical noise segments, and 10 other clustered detections were not verified as earthquakes. In total, 432 clustered events were therefore confirmed as earthquake events by visual inspection, 159 events of which (about 36%) occurred during night hours. In Fig. 6b, the red circles represent the 432 confirmed events that are distributed from the July 31 to September 18, spanning almost the entire processing time. The six template events are distinguishable among the other detections by their time and the detection values, which are 3 for all of them. We also see that the largest event in the sequence, T5, is not among the re-clustered events.
In the time period for which the BGS reports earthquakes associated with the PNR site, between August 15 and September 14, the total number of detections Some of the events missing in the BGS catalog are detected by the TM approach with detection values above 2 (maximum value of cc is 3). In addition, there are many confirmed events with low detection values. It should be noted that total number of true events are not limited to the 432 confirmed events, and among the non-clustered events there are also small events, however we could not check all of them case by case. Considering the magnitude of template events as references, the relative magnitude (Schaff and Richards 2014) of the confirmed events are estimated to be above −0.5 M L with uncertainty of ±0.5 M L . For each event, all templates are used to estimate the relative magnitudes and the uncertainty is estimated by considering standard deviation of all estimated values. Figure 7 summarizes the results of the waveformbased re-clustering of all detections. Events in clusters with more than 10 detections are shown. Cluster C3, with 115 events, is the most populated cluster. Events in this cluster happened between August 15 and September 14. Among them 38 events (33 percent) are reported events by BGS. In addition, template T6 belongs to this cluster and the majority of the events in this cluster happened before T6. Cluster C5 with 69 events is the second largest cluster and 47 of those events are also reported by BGS. Template T1, T2, T3 and T4 are in this cluster. Events of this cluster occurred between the August 21 and September 14.
The events in cluster C6 are visually verified as small seismic events, with one reported in the BGS catalog.

Detections at other stations
Waveforms of other more distant stations, such as AQ03, AQ05, AQ06, AQ07, L001, L002 and L009 have been scanned using the single station template matching algorithm. The number of detections and common detections with BGS catalog decreases dramatically for the furthest stations, such as L009 and AQ07. Detections with at least one picked P-phase are selected from all stations and are merged if their arrival time difference matches with the theoretical values. For 98 events, more than two P-phase pickings (different stations) are found. Results are shown in Fig. 8, with stations sorted based on their distance to the epicenters. We study the activated fault mechanism by waveform modeling for events whose signal quality at an appropriate number of stations is high. Source modeling for pure deviatoric moment tensors is performed using full waveform information. The inversion is based on an oc-tree grid search approach (Lomax and Curtis 2001) solving a tensor combination method (Krizova et al. 2013). Sampling is done in a 4D parameter space for strike, dip, rake and compensated linear vector dipole (CLVD) components (Lay and Wallace 1995;Aki and Richards 2002) following the uniform parameterization by (Tape and Tape 2015). The 1D velocity model for the UK (Booth et al. 2001) is used and a synthetic green's function database is computed by the Fomosto extension of the Pyrocko framework ). Our approach allows for nonlinear inversion taking into account uncertainties on the underlying source parameters.
Moment tensor solutions for the largest events were derived in two steps: (1) solving for the predominant pure double couple (DC) mechanism in the full parameter space and (2) including a CLVD while allowing for small changes of the DC part (3D+1D). During the inversion, we tolerate small truncated station-wise time-shifts to account for path effects and erroneous locations. Output results contain uncertainties of the four selected source parameters within 10% of the highest probable solution and the moment magnitude (compare Figs. 9 and 10).
Event selection is based on the initial local magnitude taken from the BGS catalogue (M L ≥ 1.0), the signal quality in the analyzed frequency band and number of usable traces. Selection of the traces was performed manually considering signal to noise ratio and quality among each component in comparison to nearby stations. We only inverted for events that show clear signals at a minimum of four stations and at least 6 traces. The BGS event locations are considered in modeling process. With an exception of T6, we were able to invert the source mechanism for all template events. Furthermore, only results with a variance reduction (VR) exceeding 50% in a 4s (T1-T3) and 30 s (T4 and T5) signal time interval are used for interpretation. Events T4 and T5 were further examined in different frequency bands and source depths and the first motion polarities were included. T4 contains no source signal below 0.3 Hz and returns unstable solutions for frequencies larger than 0.85 Hz. T5 excites signals up to 0.15 Hz while becoming unstable for frequencies larger than 0.75 Hz. Inversion results in Fig. 9 Inversion results including uncertainty information and waveform fit of selected stations for event T4. Black lines represent the observed displacement seismograms and red ones the synthetics solutions of the inverted source mechanism within 10% of the highest probability. Traces highlighted in gray were not used during the inversion We investigated, within the 4-dimensional parameter space, the best solution and, if existing, multiple minima that would indicate alternative mechanisms. Our tests show that the inversion results are stable and unique solutions of all examined events are found. The resulting 1D moment tensor solutions (Table 5 and Figs. 9 and 10) display predominantly strike-slip solutions (see Fig. 1). Observed strike directions are 15 • ± 4.5 • (T1-T4), with an exception of the largest event (T5) at around 30 • . With a dip angle of 81 • and 85 • , event T4 and T5 are steeper than T1 to T3 with dip angles between 60 • and 64 • . The rake angles for all five events are between −14 • and −7 • . Overall, events T1 to T3 exhibit large similarities in their derived source mechanism with a Kagan angle of 3.9 ± 0.8 at a variance reduction (VR) between 0.53 and 0.65. T4 and T5 reach a VR of 0.76 for a full deviatioric tensor. For the inversions of the three largest events (T1, T4 and T5) we allow for additional CLVD contributions as we deemed a sufficient enough signal quality on all three components. Smaller events (T2 and T3) show limited, or mostly low, signal quality on the vertical traces, which forced us to solely rely on horizontal recordings. As a result, we did not invert for the CLVD components of smaller events.

Event location
The aim of event location undertaken herein is to support the recognition of the active and auxiliary fault plane based on the spatial distribution of events, and to distinguish between right and left lateral mechanisms. High-quality phases obtained automatically using template matching algorithm and cross-correlation are used along with the hypoDD algorithm to carry out double-difference relocation (Waldhauser 2001). The double-difference equations are solved using singular value decomposition. For all pairs of events recorded at each station, where each pair includes a verified event and a template event, differential travel times of P-and S-phase are calculated using cross-correlation. The hypoDD method requires a weighting factor for individual differential travel time data in order to take into account the confidence in each measurement. The cross-correlation coefficient between the template's waveform and each event is used as weighting coefficient for each differential travel time. A minimum threshold value of 0.8 (out of 1) is taken for filtering out the less accurate data. In addition, dynamic weighting is necessary to exclude outlier and large (beyond 0.7 km) inter-event separations. The initial locations of all events are assumed to be at the location of the biggest event, T5 and the Preese Hall Vp model is used (Eisner et al. 2011). Waveforms of 6 stations close to the event's epicenter, i.e. stations AQ03, AQ04, AQ05, AQ06, L001 and L009 are used initially to calculate differential travel times. Stations L001 and L009 only include waveforms of the template events. Only events with at least 6 picks are included in the relative location process. P-and S-phases data are considered with equal weights.
We could calculate relative locations for 26 events using the selection and weighting criteria detailed above. The locations relative to the cluster centroid are plotted in Fig. 11. The date of occurrence of the earliest event is August 16, and for the latest is September 12. Epicentral errors in easting direction are larger than northing direction for most of the events due to  In Fig. 11b, the strike orientation of two fault plane solutions is shown, while in Fig. 11c, the depth cross-section of events is shown. Figure 11d uses beach-ball plots to represent the depth cross-section of the focal mechanisms of the 5 template events. The dip angles of the fault planes in the depth-east cross-section are shown according to each focal mechanism solution. The depth-easting pattern of the event locations can be seen to match with the calculated dipping angle of the focal mechanism solutions, differentiating the active fault plane.

Interpretation and discussion
Results of TM are compared to the BGS catalog and the catalog provided by the operator (Kettlety et al. 2020) to make conclusion about the capabilities of surface stations in shallow depth induced microseismicity studies. It is worth mentioning that the BGS locations are generally obtained on the operator surface network, with several stations within a few kilometers of the epicentre that are not considered in this study. In addition, temporal changes and characteristics of TM detections as well as faulting and source mechanism of larger events in the studied induced seismicity are discussed in this section.

Comparing TM detections with other catalog events
Assuming BGS catalog as a benchmark, results of TM detections at the site AQ04, about 1.3 km from the events' epicenters, is satisfactory as 98% of cataloged events are among the detections. In addition, the number of verified detected events using TM at AQ04 is at least 3 times larger than cataloged events for the same time span. On the other hand, the total number of detections and number of shared events with BGS catalog are rather small for the stations at greater epicentral distances. For instance, applying a bandpass filter of 20-50 Hz, and using MAD 9 as threshold value, just 12 cataloged events are detected at station L009 (Table 4). This is due to decreasing amplitude of small magnitude events with increasing epicentral distance along with high and variable noise levels.
The noise PPSD (McNamara and Buland 2004) calculated for the 3 components of all recording stations (Fig. 2) shows that station AQ04 experiences the lowest noise level in the high frequency bands, which, along with the close proximity to the event epicentres, explains the highest number of common detections with BGS. In the high frequency band, stations AQ03 and AQ05 show noise level differences, with AQ03 (nearer to the coast) slightly noisier than AQ05. However, the number of detections at AQ03 is larger than at AQ05, while the number of common events with BGS catalog is larger at AQ05 in 2-20 Hz and 20-80 Hz frequncy bands (see 5th column of Table 4). This suggests that proximity to the events and observed variability in noise level are both important factors in event detection results using the TM technique. The temporary stations of LV network show higher noise level than the other stations (which are semi-permanent). This can explain lower number of valid event detections at LV sites (see Table 4). Kettlety et al. (2020) reported in total 55555 microseismic events (mostly above M W > −2) using a dense down-hole array at about 2 km depth and an automatic event detection technique, about 135 of those events are larger than 0.0 M W and might be observed at the surface. We confirm that 96% of larger events of operator catalog are detected by TM technique. The operator catalog does not include the location of larger event T5 (2.9 M L ), since at the occurrence time of T5 (3 days after the pumping was terminated) the downhole array were not operating. From 1590 detections with TM (since August 13), 815 events could be confirmed through comparison with the operator catalog (about 50%). Their magnitude ranges between M W −2.3 and M W 1.1. The magnitudes estimated for template events (except T2) are lower in the operator catalog compared to the BGS catalog (see Table 1). This might indicate that the employed magnitude estimation method by the operator underestimates the magnitude of the most events.

Temporal changes and characteristics of TM detections
Hydraulic fracturing of the horizontal well at a depth of 2.1 km, with horizontal extension of 750 m, at this site began on August 15, 2019, and was suspended on August 23, 2019 (OGA 2019; Cuadrilla 2019). The template matching technique shows weak detections happening prior to this fluid injection and all cumulative plots (Fig. 5) show jumps on August 8, probably related to well/flow testing, followed by almost constant cumulative numbers until August 15. We could verify some of these pre-injection detections as seismic events by re-clustering and visual inspection of their waveforms (red circles in Fig. 6b). Confirmed events are distributed from the July 31 to September 18, over the entire processing time, suggesting the existence of background seismicity prior to, and throughout injection. However, lower detection values (assuming 1.6 as the cross-correlation threshold where the maximum value is 3.0) for events happening before August 17 points out that the related waveforms are not well-correlated with template events, suggesting they might be contaminated with noise and are therefore more likely to be smaller magnitude events. During fluid injection, jumps in the cumulative detection curves obtained with each template are in agreement with the injection times (see vertical gray lines in Fig. 5). This indicates that the change in detection rate at a surface station can be of physical meaning and be linked to the process that is happening in the sub-surface, down to about 2 km in this case.
After occurrence of T1, the first 'larger' event in the sequence (M L 1.6), the detection rate with T6 and T5 decreases while the number of detections with the other templates increases. This behaviour indicates a change in characteristics of the events being detected after the initiation of larger events. The change in character and size of events may, therefore, be indicative of a change in the source of seismicity.
The characteristics of templates in the time and frequency domain, in particular templates T5 (M L 2.9) and T6 (M L 1.0), explain their different temporal detection rates. T5 does not correlate with any event in the low frequency band (2-20 Hz) as its strong low frequency content (see Fig. 4) is unique and is even different from the other template events. On the contrary, in the higher frequency range, above 20 Hz, T5 does contribute to the new event detection task. Although T5 and T6 show no waveform correlation to one another in 20-80 Hz band, the cumulative detection curves obtained from those templates, during the time period when the fluid injection was activated, show similar trends (Fig. 5).

Faulting and source mechanism
The fault plane solutions calculated for 5 template events (Table 5) imply strike-slip mechanisms in agreement with studies of previous sequences near Blackpool (Baptie 2018;Clarke et al. 2014;Clarke et al. 2019). The predominantly negative CLVD of up to −21% for the largest event T5 and −17% for the second largest event T4, might suggest a closing or collapsing mechanism due to changes in fluid pressure or an additional non-strike-slip event, presumably a normal faulting minor couple, with similar strike and dip angle.
The microseismic events locations (SNR > 20), according to the operator catalog, are shown in Fig. 12a (map view) and Fig. 12b (cross-section), highlighting those with magnitude above M W 0.6 (red) along with the template events (blue). The nodal lines, calculated for each template event are shown in this figure as well. Although events location indicate a general northwest-southeast direction, the alignment of T1-T4 in a northeast-southwest direction is verified in agreement with relative location we obtained with hypoDD (Fig. 11). This direction is almost perpendicular to the horizontal wellbore profile (Cuadrilla 2019) and can indicate the strike direction and, accordingly, the left-lateral movement of the activated fault causing the larger events. The clear pattern of decreasing depth of events toward the east (Fig. 11c) can also support the northeast-southwest strike direction, as the dip angle related to that (solid lines in Fig. 11d) shows better agreement to the pattern of events. However the number of events might not be sufficient to make a certain conclusion. In the cross-section view (Fig. 12b) the pattern of larger events is also more in favour of west-dipping nodal planes, indicating the left-lateral mechanism.
T5 and T6 are located noticeably outside the main seismic cluster according to the relative location we obtained using hypoDD (Fig. 11a) and the operator's microseismic catalog (Fig. 12a). According to hypoDD, T5 and T6 locate to the west with respect to the other templates (Fig. 11a) and at a greater depth (Fig. 11c). The location of M L 2.9 event, T5, is estimated from just the surface arrays and verified based on spatial distribution of its aftershocks that occurred on a plane with northwest-southeast alignment (Kettlety et al. 2020).
At least two distinctive sequences of clustered events are observed (Fig. 7). Each of these clusters include at least one template event. Template T6 is included in the largest cluster, C3, with 115 events occurring between August 15 and September 14. Templates T1, T2,T3 and T4 fall in another sequence, C5, which includes 69 events and starts on August 20 and lasts until the September 14. The different characteristics of T6 from the other templates could be linked to the location difference, with T1-T4 belonging to a fault different to T6. Comparing the locations of events in clusters C2, C3 and C4 based on the operator's micro-seismic event catalog, we clearly see that the waveform clustering also indicates spatial clustering. Events of cluster C5 are separated from two other clusters both in map and depth view (Fig. 13). This suggests that T6 may be activated by fracture propagation, while other templates are generated by the reactivation of an existing fault, or fracture propagation in a different direction.

Conclusion and outlook
In this study, our aim was to examine the induced microseismicity related to the August 2019 hydrofracking that took place near Blackpool, UK. We used continuous seismograms recorded at a sparse network of surface broadband stations; the closest station to the fracking zone, AQ04, was located at a distance of only 1.3 km away. We concluded about the capabilities of surface stations alone to detect and characterise the induced seismicity, in terms of number of detected events, temporal variation of detections, faulting and source modeling and moment magnitude of larger events.
About the event detections and temporal variation of detections, we observed and concluded the following: (1) The number of verified events at a single surface station is at least 3 times larger than the reported events by BGS for the same time span. (2) Numerous small events prior to the main stimulation indicates the background seismicity related to the existing faults and the pre-stimulation well/flow testing. (3) During hydraulic fracturing, changes in event detection rate at the surface station correlate with the timing of pumping stages. (4) Weeks after the fracking was terminated, microseismic events continued to Fig. 13 Location of the events in cluster C2, C3 and C5, in map view (a) and depth view (b) based on the operator catalog. For more explanation please see Fig. 12 occur, and waveforms of later events show similarity to both clusters of events which were activated during the injection (cluster C3 and C5). 5) Comparing of the verified events detected by TM and the events detected and located by the operator's microseismic down-hole network, we observed that 96% of larger events (M W > 0) of operator catalog and 50% of events with magnitude above −0.5 are detected by TM technique.
By analysing the larger events and from location of the clustered events, we concluded about faulting and source mechanism of the induced seismicity. According to the waveform similarity measurement, aligned location and focal mechanism solution, T1, T2, T3 and T4 can be related to a left-lateral strike-slip fault system with northeast-southwest striking direction between 7 • and 18 • and dipping angle between 60 • and 85 • . T5 and T6 are located noticeably outside the main seismic cluster according to the relative location we obtained using hypoDD and the operator's microseismic catalog (Kettlety et al. 2020). T5, the largest event, does not fit in any TM clustering group and shows limited correlation to the remaining events and potentially activates a much larger area than T1-T4. However, all derived mechanisms, in agreement with mechanisms related to the previous sequences in 2011 and 2018, show similarities in their major double couple part. Accordingly, we can conclude the leftlateral fault mechanism for T5 as well. T1-T5 can be related to the reactivation of a same fault or different faults close to each other with similar strike direction. The predominantly negative CLVD for the largest events leads us to suggest a collapsing or even a potential minor normal faulting mechanism additional to the major strike-slip movement, which might be related to the pressure release that took place after the termination of the experiment and before the occurrence of events T4 and T5.
The moment magnitude calculated for 5 templates is less than the local magnitude estimated by BGS (see Table 1). The moment magnitude for the largest event is consistent with the M W used in the risk analysis by (Edwards et al. 2021).
The main conclusions and outlook of this paper can be summarize as: -The repeating nature and waveform similarity of fracking-related induced seismicity is perfectly suited for the application of match filter technique to detect multiple events by just using a single surface station. -Temporal behavior of earthquake driving forces can be studied by waveform-based event clustering technique. -Changes in event detection rate at a single station at the surface can give insights about the ongoing process at the subsurface giving rise to microseismicity. -In the case of occurrence of natural swarms or increasing microseismicity with unknown reasons, template matching techniques can be of great value in event detection task in absence of a well-distributed array of stations.