**2. Methodology**

Stacking-based methods have been widely used to detect and locate seismic events at di fferent scales. These methods usually consider primary phases only; thus, they can also be called partial waveform stacking [33]. This category mainly includes the di ffraction stacking (DS) and cross-correlation stacking (CCS) methods. Figure 1 shows a sketch of the principles for the two methods. The generalized equation of stacking-based seismic location methods is summarized in the following Equation (1) [13]:

$$S(\mathbf{x}, t\_0) = \sum\_{\mathbf{N}} \mathcal{C}F(t) \cdot \mathcal{M}(t) \tag{1}$$

where *<sup>S</sup>*(**<sup>x</sup>**, *t*0) is the imaging value, **x** denotes the source coordinates, *t*0 denotes the origin time, and in the waveform data, *N* is the total number of receivers, *CF*(*t*) is the characteristic function (*CF*), *t* denotes the time sample, and *M* is the imaging or migration operators.

**Figure 1.** Sketch of the principles for the DS method (**a**) and the CCS method (**b**).

Equations (2) and (3) are the formulae of the DS and CCS methods [30]:

$$S\_{DS}(\mathbf{x}, t\_0) = \sum\_{i=1}^{N} CF^i(t) \delta[t - (t\_0 + t\_{i, \mathbf{x}})] \tag{2}$$

$$S\_{\rm{CCS}}(\mathbf{x}) = \sum\_{\substack{i=1 \ i \neq 1 \atop j=i+1}}^{N} \mathcal{C}^{ij}(t) \delta \left[ t - \left( t\_{i,\mathbf{x}} - t\_{j,\mathbf{x}} \right) \right] \tag{3}$$

where *SDS*(**<sup>x</sup>**, *t*0) and *SCCS*(**x**) are the stacking values of the DS and CCS methods, *CF<sup>i</sup>*(*t*) is the *CF* of the waveform from receiver *i*, *Cij*(*t*) is the cross-correlation waveforms of the *CF*s from the pair of receivers *i* and *j*, and δ[*t* − (*t*0 + *ti*,**<sup>x</sup>**)] and δ*t* − *ti*,**x** − *tj*,**<sup>x</sup>** are the imaging operators for the DS and CCS methods, where δ is the Dirac delta function and *ti*,**<sup>x</sup>** is the theoretical travel time from receiver *i* to the source **x**. As shown in the imaging operators, DS and CCS utilize different travel time information and involve different imaging patterns, which resemble spherical surface intersection and hyperboloid intersection in the 3D scenario, respectively [30]. Consequently, the DS method generally exhibits a higher vertical imaging resolution than CCS for surface monitoring. The basic workflow of stacking-based methods includes signal pre-processing (e.g., band-pass filtering and *CF* calculation), model discretization, travel time table generation, waveform stacking, source location picking, and post-processing (e.g., enhancing the imaging resolution using a probability density distribution or Gaussian filtering).

As indicated in the above equations, DS searches both the spatial location and the origin time, while CCS decouples the origin time by making use of the cross-correlograms of the waveforms from pairwise receivers. Elimination of the origin time makes the traditional 4D source imaging problem a 3D problem for seismic location in a 3D scenario and alleviates the inherent origin time-depth trade-off in seismic source location [27,34]. Generally, the origin time is a minor parameter and can be calculated posteriorly by subtracting the theoretical travel time from the observed arrival time with respect to the same event. In this application study, we estimate the performance of stacking-based methods by only evaluating the reliability of spatial source parameters. There are numerous *CF*s being studied and utilized, including waveform envelope, semblance, kurtosis, and short-term average to long-term average ratio (STA/LTA) [28,32,34–38]. The STA/LTA is used in this study, and different lengths of short- and long-term windows are adopted for waveforms recorded at different scales (see Table 1 for more details). The above stacking-based location function can be naturally treated as a global optimization problem, in which the spatial and temporal coordinates corresponding to the maximum imaging value, max4*SDS*(**<sup>x</sup>**, *<sup>t</sup>*0)5 or max4*SCCS*(**x**)5, are generally selected as the inverted source location [39]. When multicomponent data are available, one can consider incorporating information of both compressional (P-) and shear (S-) waves to enhance the constraints. However, the complex seismic phases (e.g., uneven distributions of S-wave energy due to complex source mechanisms and the acquisition geometry) and unreliable multiple velocities can potentially deteriorate the location results by introducing additional interference instead of constructive constraints [30]. Therefore, one should be careful with multiple phases, and the waveform stacking of a single phase (i.e., P-wave) is considered here.


**Table 1.** Basic parameters of the field datasets for stacking-based location.

#### **3. Multiscale Field Data Examples**

In this section, we investigate and compare the performance of the two stacking operators based on three field datasets associated with induced seismicity at di fferent scales. Table 1 lists the basic parameters for the seismic source location of these examples. More detailed descriptions of the datasets and the associated projects can be found in the following subsections and related references.

#### *3.1. The Small-Scale Example*

Small-scale hydraulic fracturing experiments can build a bridge between experimental study and field operations, improve the understanding of the role of reservoir stimulation, and enhance the applicability of results from laboratory experiments. An in situ stimulation and circulation experiment at the Grimsel Test Site (GTS) was conducted in 2017 [40,41]. A series of small volumes of water (up to 1 m3) were injected into pre-existing faults to induce rock fracturing, which was accompanied by abundant microseismic events. A 32-channel acquisition system was used to record microseismic signals with a 1 MHz sampling rate. Here we tested the two stacking operators on field microseismic events recorded during fracturing tests in SBH-3. The homogeneous model with a P-wave velocity of 5150 m/s and band-pass filtered (1–50 kHz) waveforms from vertical components were used to locate these events. The raw waveforms of two sample events and the imaging results are shown in Figure 2. The two-dimensional (2D) imaging results are three orthogonal slices intersecting at the location of the maximum imaging value.

**Figure 2.** Waveforms and imaging results of two sample events. (**<sup>a</sup>**,**b**) Raw waveforms of the vertical components of two sample events; (**<sup>c</sup>**,**d**) 2D slices from the DS method for the events shown in (**<sup>a</sup>**,**b**); (**<sup>e</sup>**,**f**) corresponding results from the CCS method. White reversed triangles denote sensors, and white circles denote the reference locations.

Although the source energy of most events is well focused, the imaging resolutions of the two waveform stacking methods for many events are low and show obvious band-like artifacts in the east direction due to the irregular and limited coverage of the sensors. The white reversed triangles in Figures 2–4 denote the sensor locations, which mainly cover the north and depth directions and overlap in the east direction. Resultantly, the source energy is focused well in the depth–north profile, but there are strong artifacts along the east direction, indicating large uncertainty of the locations (Figure 2f) and leading to large horizontal bias of the location results when compared with those by joint hypocenter determination in a previous study [41] (Figure 3a,c and Figure 4a,c). Table 2 lists the detailed location biases of these examples. The sensor network here resembles a joint surface and downhole monitoring array, which provides good constraints for the depth and yields an average depth bias smaller than 2 m (see Figures 3b and 4b). Other factors contributing to the location uncertainty include the lack of sensor calibration and the simplified homogeneous and isotropic velocity model, which involve arrival time and travel time errors.

**Figure 3.** Location results of the DS method for the selected events. (**<sup>a</sup>**–**<sup>c</sup>**) 2D slices; (**d**) 3D view. Reversed triangles denote sensors, red dots are the reference locations, and blue dots are results from the DS method.

**Figure 4.** Location results of the CCS method for the selected events. (**<sup>a</sup>**–**<sup>c</sup>**) 2D slices; (**d**) 3D view. Reversed triangles denote sensors, red dots are the reference locations, and blue dots are results from the CCS method.


**Table 2.** Average location biases of the two stacking-based methods at multiple scales.

#### *3.2. The Exploration-Scale Example*

Seismicity related to underground mining has been observed since the beginning of the 20th century [42]. Seismic activities are closely related to the safety of mining operations, and rock bursts are often the main cause of mining accidents [43,44]. Passive seismic/microseismic monitoring can effectively detect and evaluate the seismic activity surrounding underground mines and provide early warning of potential geological risks [3]. From June 2006 to July 2007, a temporary network (HAMNET) was set up to monitor mining-induced seismicity in the Ruhr area, Germany [45]. The network consisted of 15 three-component surface stations, covering an area of about 3 km × 2 km. The depth of the longwall was about 1 km. We selected 200 events from the dataset to test the two stacking operators, and only P-waves in the vertical components were used in the location process [30]. The sampling time was 5 ms and the raw waveforms were preprocessed by a band-pass filter (5–30 Hz). We set the

target location volume as 5 km × 5 km× 5 km with grid spacing of 50 m. The imaging and location results are shown in Figures 5 and 6.

**Figure 5.** Imaging results of a sample event. (**a**) 2D slices from the DS method; (**b**) corresponding results of the CCS method. White reversed triangles denote sensors, and white circles denote the reference locations. The reference point (east, north) = (0, 0) corresponds to east = 411,117 and north = 5,721,611 in the Universal Transverse Mercator (UTM) coordinate system.

**Figure 6.** Location results of selected events. (**a**) Comparison between reference locations and results of the DS method; (**b**) comparison between reference locations and results of the CCS method. Red dots are the reference locations, and blue circles are results from stacking-based methods.

As shown in Figure 5, the horizontal resolutions of the imaging results for both methods are good due to the good horizontal coverage of the 15 surface stations. As mentioned before, the vertical imaging resolution of DS is higher than that of CCS in the case of surface monitoring. Figure 6 shows that the overall distributions of the selected microseismic events determined by the two methods are in good agreemen<sup>t</sup> with that of travel time inversion, and the CCS method has better event clustering and smaller location bias (see Table 2) than the DS method. Please note that the imaging resolution addresses the quality of the imaging profile and the location uncertainty describes the stability of the location results. The two parameters are not necessarily consistent, though they are closely related to each other. In this case, the CCS method produces more reliable and clustered location results, though it has lower imaging resolution.

#### *3.3. The Regional-Scale Example*

Induced seismicity monitoring is of grea<sup>t</sup> significance to both industrial production and public safety. During the past decade, several large-magnitude earthquakes related to hydraulic fracturing in unconventional reservoirs have been reported in the United States, United Kingdom, Canada, and China [7,9,46]. In particular, the number of seismic events observed in Oklahoma, the United States has increased significantly since 2008. Besides this, seismic activity related to the development of Mississippi limestone reservoirs in central and northern Oklahoma and southern Kansas has been occurring. In 2016, more than 1800 vertical component nodal seismometers, named the LArge-*n* Seismic Survey in Oklahoma (LASSO) array, were deployed in Grant County, Oklahoma, to study induced seismicity associated with production of the Mississippi limestone play [47]. The LASSO array covered an area of 25 km × 32 km. The sampling rate was 500 Hz and the grid spacing was set to 400 m (referring to the location error in the catalogue). We selected two events from the catalogue (event no. 23196, local magnitude (ML) 3 and event no. 23897, ML 1.5) and used a layered isotropic velocity model to image the sources [48]. The layout of the LASSO array and partial waveforms of the two sample events are shown in Figure 7. There were 1825 and 445 e ffective traces for event 23,196 (ML = 3) and event 23,897 (ML = 1.5), respectively. Our tests showed that waveforms from much fewer traces can focus the source energy. Figure 8 shows the imaging results using 10-times fewer traces.

**Figure 7.** The LArge-*n* Seismic Survey in Oklahoma (LASSO) array (**a**) and raw waveforms of the two events (**b**,**<sup>c</sup>**). The blue reversed triangles are seismometers, and the red dots are the locations of the two events from the catalogue. The reference point (east, north) = (0, 0) corresponds to east = 579,000 and north = 4,051,000 in the UTM coordinate system. (**b**,**<sup>c</sup>**) are waveforms of event 23,196 (ML = 3) and event 23,897 (ML = 1.5).

**Figure 8.** Imaging results of two sample events using 10-times fewer traces (183 and 45 traces in total for event 23,196 and event 23,897, respectively). (**<sup>a</sup>**,**<sup>c</sup>**) 2D slices from the DS method for the events shown in Figure 7b,c; (**b**,**d**) corresponding results from the CCS method. White circles denote the reference locations from the catalogue.

The horizontal imaging resolutions of the two methods are very high, while the vertical resolutions are much lower. The inverted depths from the two methods are quite consistent. The location bias between the stacking-based methods and the catalogue is within 2 km (see Table 2), most of which comes from uncertainty of the depth, and is acceptable considering the large scale of the monitoring array and the target volume. Although most seismic events were recorded by more than 100 or even 1000 seismometers in the LASSO array, the imaging results in Figure 8c,d indicate that several dozens of traces are sufficient to focus the source energy of these induced earthquakes, as long as a good spatial coverage of the selected traces can be ensured.

#### *3.4. Comparison of the Results Across Scales*

Compared with the DS method, the CCS method makes use of more waveform redundancy and includes more constraints for the source imaging process, which is beneficial for enhancing the signals and suppressing the artifacts. However, CCS also introduces more potential noise energy and interference information, which naturally reduces the imaging resolution (see Figures 2, 5 and 8). The results also demonstrate that the imaging resolution is highly dependent on the acquisition geometry. For the small-scale example, the sensors surround the stimulation site and ensure a relatively good imaging resolution for both the horizontal and depth directions (Figure 2c–e), though partial events have a large horizontal bias due to the limited coverage of sensors in the east direction (Figure 2f). For the exploration-scale example, the horizontal coverage of the network is two-times the target depth, which can produce reliable depth locations (Figure 6) [49]. For the regional-scale example, the vertical

imaging resolutions of the two sample events are poor, due to the events being located near the margin of the array (Figures 7 and 8).

Table 2 shows the average location bias for the three examples, which was calculated by dividing the total absolute bias between the location results of stacking-based methods and the reference locations by the number of selected events. Figure 9 shows histograms of the overall location biases for all selected events in the small-scale and exploration-scale examples. The comparison of the results across the scales clearly shows that the location bias is consistent with the scale size, which naturally results from the different frequency contents of the seismic waveforms and grid spacing values (see Table 1) at different scales.

**Figure 9.** Histograms of the overall location biases of the DS and CCS methods for the small-scale (**a**) and exploration-scale (**b**) examples.

#### **4. Discussion and Conclusions**

In this work, we investigated the feasibility and potential of stacking-based methods using production-induced seismicity at different scales. The examples revealed the characteristics of the imaging resolution for the two stacking operators and demonstrated their feasibility in locating seismic events at multiple scales. A comparison of the location results across the scales indicated that the imaging resolution is highly dependent on the acquisition geometry, and the location bias is closely related with the scale. Although the reference locations determined by other methods or from the catalogue may not be reliable, the comparison can still indicate several potential factors that are responsible for the accuracy and reliability of the methods—that is, the layout of the monitoring network, the reliability of the velocity model, and the quality of the waveforms. The methods are noise-resistant and automatic, which means that weak events with low signal-to-noise ratio can be located and no phase picking procedure is needed.

So far, field applications of stacking-based methods have mainly focused on the exploration scale—e.g., for hydraulic fracturing monitoring with large and dense arrays in unconventional oil and gas reservoirs [27,34]. More recently, these methods have been naturally adapted to analyze low-frequency seismic events at larger scales [28,50], where a large depth uncertainty exists. There are also emerging applications for seismic events with comparably high-frequency content at smaller scales—e.g., rock burst and acoustic emission, where the issues of location uncertainty and velocity reliability are matters of considerable concern [32,51]. The events tested in this study have relatively high signal-to-noise ratio, and waveforms recorded with dozens of receivers can recover the source energy quite well. Further work will consider weaker seismic events associated with hydraulic fracturing in unconventional and geothermal reservoirs and foreshock/aftershock activity potentially preceding/following tectonic earthquakes. It is worth noting that stacking-based methods are more resistant to white noise, and one should use an advanced filtering technique to address the influence of correlated noise from, for example, vicinity to hydraulic fracturing wells [51]. For events with small magnitudes and/or low signal-to-noise ratios, wavefield reconstruction and enhancement may be required to ensure a reliable source location when using waveform stacking techniques, especially when partial traces are noisy or disabled. Another challenge for stacking-based location methods is the relatively high computation cost due to waveform storage and transmission, especially for real-time applications with large monitoring arrays. High-performance computing facilities [52,53] and advanced inversion algorithms [35,39] are promising solutions for improving the computational e fficiency.

The current study is intended to be a starter for more in-depth and comprehensive investigations of stacking-based methods across scales, and the conclusions from this study still remain at a relatively qualitative level. Specifically, the stacking operators should be studied in a more systematic manner by adjusting the characteristic functions and pre-filtering techniques of the waveforms, acquisition geometries/channels, and velocity models. Although there have been several attempts at comparing the performance of di fferent characteristic functions [17,18], field applications at multiple scales are more data- and problem-dependent. Meanwhile, novel and potential stacking operators producing better seismic energy focusing are also worth investigating. Further benchmarking studies are required to improve the applicability of stacking-based methods and build a suitable workflow for seismic location at multiple scales.

**Author Contributions:** Conceptualization, L.L., Y.X. and J.T.; methodology, L.L. and Y.X.; software, L.L.; validation, L.L.; writing—original draft preparation, L.L.; writing—review and editing, Y.X. and J.T.; supervision, J.T.; funding acquisition, L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Natural Science Foundation of Hunan Province, gran<sup>t</sup> number 2019JJ50762, the Open Research Fund Program of State Key Laboratory of Acoustics, Chinese Academy of Sciences, gran<sup>t</sup> number SKLA201911, and China Postdoctoral Science Foundation, gran<sup>t</sup> number 2019M652803. The work of Y.X. is funded by the Natural Environment Research Council, gran<sup>t</sup> numbers NE/M003507/1 and NE/K010654/1, and the European Research Council, gran<sup>t</sup> number GA 638665.

**Acknowledgments:** We acknowledge the open source datasets from the GTS experiment, the HAMNET network, and the LASSO array.

**Conflicts of Interest:** The authors declare no conflict of interest.
