Application of Waveform Stacking Methods for Seismic Location at Multiple Scales

: Seismic source location speciﬁes the spatial and temporal coordinates of seismic sources and lays the foundation for advanced seismic monitoring at all scales. In this work, we ﬁrstly introduce the principles of di ﬀ raction stacking (DS) and cross-correlation stacking (CCS) for seismic location. The DS method utilizes the travel time from the source to receivers, while the CCS method considers the di ﬀ erential travel time from pairwise receivers to the source. Then, applications with three ﬁeld datasets ranging from small-scale microseismicity to regional-scale induced seismicity are presented to investigate the feasibility, imaging resolution, and location reliability of the two stacking operators. Both of the two methods can focus the source energy by stacking the waveforms of the selected events. Multiscale examples demonstrate that the imaging resolution is not only determined by the inherent property of the stacking operator but also highly dependent on the acquisition geometry. By comparing to location results from other methods, we show that the location bias is consistent with the scale size, as well as the frequency contents of the seismograms and grid spacing values.


Introduction
Seismicity can occur naturally in seismogenic areas or be induced by industrial operations, ranging from laboratory acoustic emission events to large-scale global earthquakes. Under controlled conditions, laboratory and small-scale experiments can reveal the mechanism of fracture initiation and propagation, quantify the changes of reservoir permeability, and evaluate the stiffness of natural fractures or faults [1,2]. At the exploration scale, passive seismic monitoring has been used to delineate fracture propagation, monitor reservoir deformation and fluid migration, and assess seismic risks associated with many subsurface operations, such as underground coal mining, hydraulic stimulation of unconventional oil and gas reservoirs, geothermal exploitation, and carbon dioxide injection and geological storage [3][4][5][6][7][8][9]. A complete and accurate earthquake catalogue is an important basis for subsequent data processing and even directly determines the reliability of seismic monitoring and earthquake prediction at the regional or global scale [10]. Seismic source locations are key information of earthquakes and play an important role in characterizing the geometries of multiscale fractures/faults, evaluating seismic activities, and inverting the source mechanism and in situ stress state. For instance, the spatial and temporal distribution of seismic events can help reveal the mechanism and propagation of rock fractures at the laboratory/small scale, as well as provide important information for the assessment of tectonic and volcanic seismicity at local and regional scales. Seismic location as a typical inverse problem, covering geophysical, seismological, acoustic, and engineering applications at multiple scales, has experienced significant methodological and application progress during the past century [11][12][13][14][15].
With the development of modern seismic instrumentation for dense acquisition and induced (micro) seismicity monitoring technology, new challenges and opportunities have emerged for noise-resistant, automatic, and real-time seismic location methods. Specifically, a new category of waveform-based location methods (e.g., waveform stacking methods and time reverse imaging) has emerged as a counterpart of conventional travel time-based inversion [13,[16][17][18]. Waveform-based methods locate the seismic source by combing the travel time, amplitude, and phase information of seismic waveforms or wavefields to reconstruct and focus the source energy into an image profile [13]. There are three important advantages for waveform-based seismic location methods. First, they are noise-resistant, since multichannel waveforms are involved, and the coherence is enhanced to detect and locate more events. Second, the methods are basically automatic and data-driven, which enables a more efficient location process and avoids potential subjective interference, such as phase picking. Third, the source locations are resolved as images instead of simple dots, offering more insights into source processes and surrounding structures. Stacking-based methods are the most mature and successful methods considering their wide applications across scales. The basic principle of stacking-based location methods is reconstructing and focusing the radiated seismic energy from the source with a certain stacking operator, for example, the diffraction stacking (DS) operator [19][20][21] or the cross-correlation stacking (CCS) operator [22][23][24]. The origin of stacking-based seismic location methods can be traced back to the 1990s. Kiselevitch et al. (1991) proposed to utilize the maxima of the semblance over space, time, and channels to detect and locate microseismic events, and named the method "seismic emission tomography" [25]. Later on, passive seismic emission tomography was developed to locate microseismic events using recorded waveforms from surface arrays [26]. These methods share the same principle with DS in reflection seismology, that is, stacking the waveforms of individual receivers with the corresponding one-way travel time moveout to enhance the diffracted/scattered seismic energy. CCS is another well-established stacking-based method for source location. CCS exploits correlation waveforms corresponding to differential travel times at pairwise receivers from common events.
Stacking-based methods have been applied to field data at multiple scales, including experimental microseismic/acoustic emission events, mining-induced seismicity, hydraulic-fracturing-induced seismicity, volcanic-tectonic seismicity, and regular earthquakes [27][28][29][30][31][32]. However, more systematic benchmarking studies of these approaches across scales are still needed. Specifically, the imaging resolution characteristics and location reliability resulting from different frequency contents at different scales need to be studied further. Different from the comprehensive literature review of the methodological and application progress of waveform-based methods in Li et al. (2020) [13], the presented work is a case study and is the first attempt at investigating the performance of two specific stacking-based methods with common field datasets across scales. We aim to investigate the feasibility and reliability of the methods through mutually comparing the imaging results and validating the location results with reference locations, which can provide guidance for further evaluations and improvements of the methods.

Methodology
Stacking-based methods have been widely used to detect and locate seismic events at different scales. These methods usually consider primary phases only; thus, they can also be called partial waveform stacking [33]. This category mainly includes the diffraction 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]: where S(x, 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.  Equations (2) and (3) are the formulae of the DS and CCS methods [30]:  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 tradeoff 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 CFs being studied and utilized, including waveform envelope, semblance, kurtosis, and short-term average to long-term average ratio (STA/LTA) [28,32,[34][35][36][37][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, x , 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 Equations (2) and (3) are the formulae of the DS and CCS methods [30]: where S DS (x, t 0 ) and S CCS (x) are the stacking values of the DS and CCS methods, CF i (t) is the CF of the waveform from receiver i, C ij (t) is the cross-correlation waveforms of the CFs from the pair of receivers i and j, and δ[t − (t 0 + t i,x )] and δ t − t i,x − t j,x are the imaging operators for the DS and CCS methods, where δ is the Dirac delta function and t i,x 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 CFs being studied and utilized, including waveform envelope, semblance, kurtosis, and short-term average to long-term average ratio (STA/LTA) [28,32,[34][35][36][37][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, max S DS (x, t 0 ) or max S CCS (x) , are generally selected as the inverted source location [39]. When multicomponent data are available, one can consider incorporating information Energies 2020, 13, 4729 4 of 15 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.

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 different 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.

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 m 3 ) 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. (e,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,  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.
Energies 2020, 13, x FOR PEER REVIEW 6 of 14 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] (Figures 3a,c and 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.

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 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 . 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.  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 agreement 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.

The Regional-Scale Example
Induced seismicity monitoring is of great 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 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 Energies 2020, 13, 4729 9 of 15 good agreement 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.

The Regional-Scale Example
Induced seismicity monitoring is of great 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 effective 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.     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.

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.

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 lowfrequency 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

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 efficiency.
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 different 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.