In addition to earthquakes, a broad range of slip events and surficial processes can generate seismic signals [e.g., Nettles and Ekström, 2010;Arrowsmith et al., 2010;Obara and Kato, 2016;Poli and Shapiro, 2022]. Transient surficial processes, when effectively coupled to the solid Earth, may produce coherent seismic signals recorded over a large region [Gerstoft et al., 2008;Schmandt et al., 2013;Nishida and Takagi, 2016]. For example, large calving events in Greenland can have equivalent magnitudes of ∼M5 earthquakes [Ekström et al., 2003;Tsai and Ekström, 2007], and small landslides can produce coherent surface wave packets up to thousands of kilometers away [Okuwaki et al., 2021;Luo et al., 2023].
The high temporal resolution of such signals can offer unique insights into the dynamics of their sources [Brodsky et al., 2003;Ekström and Stark, 2013;Iverson et al., 2015;Schmandt et al., 2017;Olsen et al., 2021]. However, these exotic seismic sources often do not produce much high-frequency radiation, and standard catalogs often miss these sources because conventional methods rely on high-frequency body waves for identifying earthquakes [Bahavar et al., 2019].
For such sources, surface waves of intermediate-to long-period often have greater signal to noise ratios than body wave phases [Ekström, 2006]. However, the emergent onsets and broad waveform packets of surface waves make it challenging to accurately determine the arrival times needed for routine detection and location operation [McGuire, 2008]. Additionally, strong crustal heterogeneities modify the arrival times in complex ways [Laske et al., 2013], also hindering their use in locating seismic sources with standard methods. Array analyses of surface waves exploit waveform coherence and have proven effective in detecting and locating both earthquakes and unconventional seismic sources [Shearer, 1994;Ekström, 2006;Gibbons and Ringdal, 2006;Fan et al., 2018;Arrowsmith et al., 2022], and have led to discoveries of new classes of seismic sources not included in regular catalogs [Chen et al., 2011;Gualtieri and Ekström, 2018;Fan et al., 2020Fan et al., , 2022]], including glacial quakes, landslides, volcanic eruptions, and stormquakes from ocean-solid Earth interaction [Ekström et al., 2003;Ekström and Stark, 2013;Poli and Shapiro, 2022;Fan et al., 2019].
Previous methods have primarily focused on using Rayleigh waves, partly because they can be recorded in the vertical component, and vertical noise is almost always much less than horizontal noise [e.g., Berger et al., 2004]. Furthermore, a vertical seismometer can be considered omnidirectional, whereas maximizing a Love-wave signal necessitates determining the direction of arrival. Therefore, Love waves have rarely been used in identifying seismic sources [e.g., Nishida and Shiomi, 2012;Matsuzawa et al., 2012], although they have -2-Confidential manuscript submitted to Geophysical Journal International a unique sensitivity to horizontal forces in exotic seismic sources, and they may be less impacted by microseisms [Cessaro, 1994;Friedrich et al., 1998;Berger et al., 2004;Nishida et al., 2008].
Here we combine Love and Rayleigh waves measured with globally distributed arrays to detect and locate seismic sources. The surface waves are filtered in the intermediateperiod band, 20-50 s, and the arrays can be spatially disconnected. Our method is built upon the AELUMA (Automated Event Location Using a Mesh of Arrays) method [de Groot-Hedlin and Hedlin, 2015;Fan et al., 2018], which explores the local coherence of surface waves between adjacent stations. Our method has the advantage of not requiring a velocity model or knowing the event type. We demonstrate our approach using synthetic data and real observations, and we show that combining both Love and Rayleigh waves can lead to more detections and better locations. The new approach improves the detection of remote events at mid-ocean ridges and Greenland. The method is one of the first attempts to exploit Love waves for identifying seismic sources, and it suggests the potential future use of records from instruments with solely horizontal components, such as tiltmeters.
We describe our procedure in detail in this section before proceeding to test its performance. The AELUMA method was first introduced to detect and locate infrasound sources [de Groot-Hedlin andHedlin, 2014, 2015]. The method was later updated to locate seismic sources using body waves [de Groot-Hedlin and Hedlin, 2018] and Rayleigh waves [Fan et al., 2018]. Although the targeted signals differ, the general principle of the approach is to use the local coherence of the recorded waves to detect signals and then aggregate the detections to locate sources. In this study, we reconfigure the method by incorporating records of horizontal components. We highlight new improvements in this study and reference some methodological details to Fan et al. [2018].
With the available stations, we first construct non-overlapping three-station subarrays (triad; e.g., Figure 1) using Delaunay triangulation [Lee and Schachter, 1980]. The ideal tree-station subarray would form an equilateral triangle to optimize resolving the propagation direction. Large or small interior angles can lead to large uncertainties in the resolved propagation directions. In practice, the triad interior angles are required to be in the range of 30 • to 120 • to robustly resolve surface wave propagation directions. To maintain the local coherence, the triad side lengths are required to be in between 10 to 600 km. Triads that do not meet the requirements are not used for further analyses. This discretization is performed for daily records to adapt to changing array configurations. As we show below, this discretization strategy can flexibly combine spatially disconnected arrays across multiple continents to optimize the azimuthal coverage for remote events.
We use -p array processing to detect coherent signals crossing each triad. We first divide the continuous data into sequential 360 s windows with a 180 s increment, and then apply the -p analysis to the windowed waveforms of each triad independently. For the vertical records, the detection procedure is similar to that outlined in Fan et al. [2018]. We pairwise cross-correlate the records and obtain three differential arrival times when the signals are sufficiently coherent (with an average cross-correlation coefficient greater than 0.5). While the frequency band of interest is much lower than the Nyquist frequency, the high sampling rate (1 Hz) provides a high temporal resolution of the differential arrival times between stations. We use these differential arrival times to invert for a horizontal slowness (local phase
where the separation distance between the th and th stations isand the differential arrival time is v , , leading to a resolved wave propagation direction
With the obtained slowness, we then slant-stack the waveforms (e.g., waveforms in the gray boxes in Figure 1B). The peak amplitude of the stack is taken as the beampower ( ) and its associated arrival time is taken as the centroid arrival time ( ) of the detected signal.
Unlike the vertical component, records of the horizontal components are often incoherent across a subarray because the propagation direction might not align with the recording directions, leading to lower average cross-correlation coefficients. To improve the coherence, we rotate the horizontal records before applying the same detection procedure as described above (Equation 1). We first assume a recording direction for a triad, and rotate the two horizontal components to direction for each station (e.g., Figure 1). The rotated records are then pairwise cross-correlated, and a wave propagation direction is estimated using the same detection procedure as for the vertical component. We grid-search from 0 • to 90 • with an increment step of ( = 1 • in this study) for the rotation. The direction that maximizes the average cross-correlation coefficient is labeled as direction 1 . For each , the horizontal components are also rotated to its corresponding perpendicular direction,
to compute a second direction 2 that maximizes the average crosscorrelation coefficient in the direction range of 90 • to 180 • . Such a searching procedure is applied to each triad at every time window. The searching procedure is independent for each time window and triad, and the detection step can be parallelized to improve the computational efficiency. For example, a 16-core 2020 MacPro takes less than five minutes to finish daily detections of ∼1000 triads. The two rotation directions ( 1 and 2 ) do not need to be perpendicular to each other. After identifying the two directions ( 1 and 2 ), the associated wave propagation directions are solved as 1 and 2 (e.g., Figure 1). The two rotation directions are different from the estimated propagation directions (e.g., Figure 1). We only gridsearch rotation directions from 0 • to 180 • as their reverse directions (180 • to 360 • ) would lead to the same measurements (propagation direction and differential arrival time).
The rotation procedure is applied to each triad independently, and the resolved 1 and 2 can be different from subarray to subarray. The example triad in Figure 1A has a 1 of 16 • and a 2 of 110 • for the data shown in Figure 1B (as described below, these are synthetic seismograms), which are perpendicular (transverse) and parallel to the propagation direction (local great circle path). As shown in Figure 1B, the rotated horizontal waveforms are crosscorrelated to obtain differential arrival times, which are then used to resolve propagation directions (Equation 1). The resolved propagation directions from the horizontal components in Figure 1 are about the same as using the vertical component, confirming the feasibility of the proposed approach. Although the example triad shows that 1 and 2 are the transverse and radial directions (Figure 1), the searching scheme does not impose constraints on the wave types as the aim is to identify a rotation direction that maximizes the wave coherence.
As we show in Section 4, mixed surface waves may still have high coherence and can provide useful measurements in surface wave propagation direction and differential arrival time. If a rotation direction is within 15 • of the propagation direction ( 1 -1 or 2 -2 , difference ranging from 345 • -15 • and 165 • -195 • ), we define the waves passing a triad as Rayleigh waves, and Love waves correspond to rotation directions within 15 • of the transverse direction ( 1 -1 or 2 -2 , difference ranging from 75 • -105 • and 255 • -285 • ). The remaining waves are categorized as mixed waves, which are the superposition of rotated Rayleigh and Love waves. We can combine the two sets of horizontal measurements to enhance the detection sensitivity. -4-Confidential manuscript submitted to Geophysical Journal International
We process the horizontal and vertical components independently to obtain three sets of measurements (centroid arrival times, direction of propagations, local phase velocities) for a triad at a given time window. The processing and quality control parameters are similar to those used in Fan et al. [2022]. Ideally, the sum of the delay times, sum = 1,2 + 2,3 + 3,1 , would be zero if the signals traveled along the assumed great circle paths between the station pairs. However, three-dimensional (3D) velocity heterogeneities may cause distortions and noise may impact the consistency [e.g., Foster et al., 2014;Fan et al., 2018]. In practice, a coherent signal is detected if the sum of the delay times, sum , is less than 60 s [Cansi, 1995].
The 60 s threshold is empirically selected as the period band of interest is set as 20-50 s.
The sliding window procedure could lead to multiple detections for the same coherent signal at a triad. We also remove duplicates at each triad for each recording component if multiple detections are identified within 180 s, and we only keep the detection with the highest crosscorrelation coefficient. We then inspect the initial detections, and discard measurements if the average cross-correlation coefficient is less than 0.6, the peak power amplitude of the stacked waveforms is less than 100 , or the local surface wave phase velocity is less than
As we are interested in surface waves, we discard detections with local phase velocities outside the range of 2.5 to 5 km/s. These local phase velocity limits are selected based on tomographic models of the western and central US [Ekström, 2014].
It is worth noting that these empirical thresholds are selected for optimizing the detection performance using the US seismic data, and the thresholds should be adjusted accordingly for other regional arrays [e.g., Okuwaki et al., 2021].
Since we want to use this method to identify events missing from standard catalogs [Fan et al., 2018], we first identify measurements that are associated with cataloged earthquakes, and then use the remaining measurements from the detection step to locate previously missed sources. It is worth noting that our method may miss cataloged earthquakes, as we use coherent surface wavefields recorded by continental-scale arrays instead of highfrequency body-wave phase picks. But we do not aim to reproduce the existing catalogs, rather we hope to provide complementary identifications of exotic seismic sources or earthquakes with abnormal seismic radiation.
We first compare the detection results for each triad with predictions from globally distributed M≥2.5 earthquakes reported by the International Seismological Centre (ISC)
bulletin [International Seismological Centre, 2013]. The ISC reports several types of earthquake magnitudes, including the body wave magnitude , local magnitude , and moment magnitude . When is available, we use it as the preferred magnitude for larger earthquakes [International Seismological Centre, 2013]. Assuming an average velocity of 3.5 km/s, we use the earthquake origin times and locations to forward-compute the predicted propagation directions and centroid arrival times at each triad. The predictions are compared with the detections, and the detections are considered to be earthquake measurements if there are more than 100 prediction-detection pairs having measurements within 180 s in time and within 20 • in propagation direction. We can combine triad results for Rayleigh, Love, and mixed waves. The differences between the observed and predicted propagation directions are taken as empirical arrival angle anomalies for the source-receiver pairs. Current surface wave velocity models cannot fully explain such observed arrival angle anomalies [e.g., Larson and Ekström, 2002;Foster et al., 2014;Pedersen et al., 2015], and we record the arrival angle anomalies for M≥3.5 earthquakes and later use them to empirically correct the measurements by removing the anomalies from the misfit residuals to improve location accuracy [Fan et al., 2018]. The detections associated with ISC earthquakes are removed from further analysis.
Before determining locations, the remaining measurements are divided into clusters using a set of potential source locations. To construct potential source locations, we first create spatial model points by combining evenly distributed points across the globe (denoted as -5-Confidential manuscript submitted to Geophysical Journal International GP) and known seismic event locations (denoted as EP). The seismic event locations (EP) are drawn from epicenters of M≥3 earthquakes reported by the ISC [International Seismological Centre, 2013] from 1995 to 2020, landslides reported by NASA [Kirschbaum et al., 2015], exotic seismic sources cataloged by Bahavar et al. [2019], and glacial quakes in Greenland [Tsai and Ekström, 2007;Veitch and Nettles, 2012;Olsen and Nettles, 2017].
We begin by generating a set of sparse points by combining evenly-spaced points (GP) that are separated by 4 • (SGP) and the event points (EP). The combined points can be treated as a two-dimensional (2D) point cloud and downsampled to obtain the final sparse points (SP). To generate the SP, the 2D surface is divided into 2 • boxes, and the points (SGP+EP)
within each box are averaged to obtain a reduced number of points, which improves computational efficiency (e.g., MATLAB function, pcdownsample, as shown in Figure 2A). Similarly, we can create a set of dense points (DP) by using globally evenly-spaced points (GP) with a 0.25 • separation (DGP), and the DP points are obtained from further averaging the combined points (DGP+EP) over 0.2 • distance. Here the potential source locations are not evenly distributed (Figure 2), which scheme prioritizes regions with active seismic events (earthquakes and non-earthquakes).
For each of the SP points, we compute predicted propagation directions and record the differences between the predictions and the measurements as for the detections. We also compute a time difference as = /3.5-, where is the point-triad geodesic-minimal-arc distances. If there are over 100 measurements within 25 • of the predicted propagation directions, we associate these measurements into clusters based on their Euclidean distances of ( , ), where denotes the triad index and the Euclidean distance limit is set as 20 for clustering. This distance limit is another empirical parameter that users can test based on the array configuration. Nearby SP points may lead to multiple clusters sharing the same measurements, and clusters are taken as duplicates if more than 60% of the detections are shared.
For such cases, only the cluster with the most detections is saved for the duplicates, while the others are discarded to save computation time. For each cluster, we fit linearly interpolated 2D planes to the measured propagation directions and arrival times, respectively, and remove detection outliers that deviate away from the planes by 10 • for propagation direction or 60 s for arrival time. Additionally, the clusters are ranked by the number of detections, and clusters with more detections are prioritized to use the detections if they are shared by more than one cluster.
The clusters are sequentially examined to locate seismic sources, and we resolve the origin times after determining the locations. The location procedure first searches the sparse points (SP, Figure 2A) to determine a preliminary location and then searches dense points (DP, e.g., Figure 2B) within 4 • of the preliminary location for a final location. The location procedure is similar to that detailed in Fan et al. [2018], which aims to minimize the differences between predicted and observed propagation directions. We use the measured arrival angle anomalies as the empirical corrections of the misfit residual for the searching pointtriad pairs when the grids are within 100 km of the closest M≥3.5 earthquakes that were recorded during the observational period [details in Fan et al., 2018]. Once a location is determined, we estimate the origin time of the event and an average propagation velocity by minimizing ℓ 1 misfit between the observed and predicted centroid arrival times from the location. The associated detections of the cluster are removed from the detection pool and they cannot be used by other clusters. As described below the location procedure can be applied to detections from the vertical component and the rotated horizontal components (e.g., ,
where is the triad index, is the beampower measurement in , and is the sourcetriad distance in degree. Here the beampower is the waveform stack amplitude at , and it is different from the conventional surface wave amplitude that is used for estimating earthquake surface wave magnitude [e.g., Shearer, 2019]. With the collection of the M SE estimates for a source, we take the median of the estimates as the source magnitude estimate,
and use its standard deviation, Var{M SE }, to infer the estimate uncertainty. We describe further details in the Results section (Section 4).
We have used two datasets to demonstrate the feasibility of using Love waves to detect and locate seismic sources and to test various aspects of the methods described above.
Both of these datasets are derived from the same source, the July 6, 2019 7.1 Ridgecrest earthquake in southern California [Ross et al., 2019]. Our first dataset is synthetic seismograms calculated from the Global Centroid Moment Tensor (GCMT) solution [Ekström et al., 2012] for this event (herein denoted as RS). We compute three-component synthetic seismograms (vertical, north-south/east-west, and radial/transverse) with a 1 Hz sampling rate. The synthetic waveforms are calculated using the Instaseis method [Driel et al., 2015].
This method generates synthetic seismograms efficiently using pre-computed Green's function databases (e.g., Figure 1B). In this study, the Green's functions are from the anisotropic version of the PREM model up to 5 s [Dziewonski and Anderson, 1981] including the Alaska leg of the USArray Transportable Array, networks of the contiguous US, and European networks (see Supplementary Materials). Data from 1362 unique stations were downloaded and used this study (e.g., Figure 3). We obtained daily, continuous, three-component long-period data (sampled at 1 Hz) for each station (see Data Availability for details). The instrument responses are not deconvolved from the records as their effects are marginal for determining the locations of seismic sources [Fan et al., 2018]. Instead, we correct the record amplitudes by normalizing the waveforms using the instrument gains. The corrected amplitudes (in ) are used for estimating an apparent surface wave magnitude (M se ) after a source is successfully located. As with RS, the real data are band-pass filtered at 20 to 50 s period (0.02-0.05 Hz). We have experimented with a range of period bands, and the 20-50 s period band seems optimal for our purpose. This frequency band has a low noise level, allowing moderate size earthquakes to be well recorded globally [Shearer, 1994;Ekström, 2006;McGuire, 2008]. Love and Rayleigh waves would provide similar measurements for the surface wave propagation direction (Figure 1), and the measurements can be used uniformly in determining the event locations. For demonstration purposes, we first separately examine the location results from 1 and 2 directions for the synthetic case (RS, Figure 4A,B). However, 1 and 2 are obtained from limited ranges of searching directions, and detections from 1 and 2 are combined in practice for better locating seismic sources (RS in Figures 4C and 5B,D). Moreover, all observations, including detections from both the vertical and horizontal components, can be combined to identify seismic events (Figure 6E). In this case, all detections ( v , 1 , and 2 ) are pooled together for clustering, and the detection with the highest average crosscorrelation coefficient of a given triad at a time window would be chosen for populating a cluster. Combining detections from different components offers multiple measurements at the same triad, which allows a better azimuthal coverage of the detections and helps to provide reliable measurements when some components are noisy (Figure 4 By applying the method to the July 2019 data (the third dataset), and combining detections from all three components, we detect 336 ISC earthquakes globally; 213 of these are within the adjacent regions of the USArray and European array (region of interest, within 90 • of 65 • /-75 • in latitude/longitude; Figure 7). Each earthquake was detected by using more than 100 triads, and lowering the detection threshold would lead to identifying more events.
Here we also define robust sources as events detected and located by more than 100 triads with maximum spatial uncertainties less than 5 • (see the following section for uncertainty analysis). There are additional 128 robust sources identified within the adjacent regions of the USArray and European array using all-three components (region of interest, Figure 8).
These new seismic sources were not reported in the ISC catalog, and they may have been missed earthquakes in remote regions (e.g., Figure 9), landslides and submarine landslides (e.g., Figure 10), and glacial quakes (e.g., Figure 11).
In Figure 3D, the synthetic example of the 2019 Mw 7.1 Ridgecrest earthquake (RS)
shows that the M SE estimate is 6.5 ± 0.2, 6.4 ± 0.2, and 7.5 ± 0.3 for using the vertical, radial, and transverse components, respectively. The vertical and radial components lead to an underestimated magnitude comparing to the moment magnitude, and the transverse component leads to an overestimation. For simplicity and consistency, we opt to use the vertical component beampowers (or the vertical measurements in the all-three components set) to estimate the detected source magnitudes (Figure 5C). For example, the 2019 Mw 7.1 Ridgecrest earthquake is estimated to have a M SE of 6.8 ± 0.3 using real data. To scale the M SE estimates to moment magnitudes, we use 71 earthquakes occurred in July 2019 that have moment magnitudes ( ) and were also detected by combining measurements from all three components (Figure 7C,F). These earthquakes are in the moment magnitude range of 3.9 to 7.1 (Figure 12). We identify two coefficients of a linear scaling relation (Figure 12B) between the moment magnitude and the M SE estimate by minimizing the ℓ 1 misfit and obtain
We use the measured vertical beampowers and the scaling relation to estimate the moment magnitudes for the newly identified events in this study (Figure 8).
To obtain coherent measurements using the horizontal components, the proposed rotation scheme is critical. Simply using the north-south or east-west components would miss many events (e.g., Figure 6). For example, we can detect 18 ISC earthquakes and 5 new sources on July 6, 2019 using the vertical component (Figure 6A), 10 ISC earthquakes and 3 additional sources using the north-south or east-west component (Figure 6B,C), 20 ISC earthquakes and 9 new sources using the horizontal components (Figure 6D), and 32 ISC -8-Confidential manuscript submitted to Geophysical Journal International earthquakes and 19 new sources using detections from all-three components (Figure 6E).
Specifically, using the rotated horizontal components, we can detect and locate seismic sources that did not generate coherent surface wavefields in the vertical component (Figure 6D). For example, an 5.4 aftershock occurred shortly after the 2019 7.1 and an 5.1 Aleutian earthquake were missed by the vertical detections but can be identified using the horizontal records (Figure 6F,G). Combining all-three components can lead to locations that cannot be resolved using only the vertical or horizontal components (Figure 6E).
In the region of interest (within 90 • of 65 • /-75 • in latitude/longitude; Figure 7), solely using Rayleigh waves (vertical component) would lead to identifying 150 ISC earthquakes during July 2019 (Figure 7A,D), and using the horizontal components (Love and Rayleigh)
waves would lead to identifying 142 ISC earthquakes (Figure 7B,E). During the same observational period, we can identify 213 ISC earthquakes by combining measurements from all-three components (Figure 7C,F), showing a significant improvement in the total detection number. Statistically, the detected earthquakes have a magnitude of completeness ( ) around 4.5 for using measurements from all-three components. The is estimated using the maximum curvature method [Woessner and Wiemer, 2005], and the estimate is the same for using the vertical, horizontal, or all-three components, respectively (Figure 7). The associated -values are obtained using their maximum likelihood estimates [Aki, 1965] with values as 1.09, 1.01, and 1.15 for using the vertical, horizontal, and all-three components, respectively (Figure 7). The different -values and the same values suggest that including all-three components may help detect more larger events, but does not significantly impact detecting smaller events. We rely on the local coherence of surface wavefields to detect and locate seismic sources, and the success of source identifications does not directly correlate with their magnitude (amplitude) but the waveform similarities.
As discussed in Section 2.2, our method takes advantage of the ISC catalog. We first use the ISC catalog to identify earthquake detections and then isolate the remaining measurements for locating new surface wave sources, which can complement the ISC catalog that is obtained using high-frequency body-wave picks. When only using the vertical component, we can obtain 69 new sources for July 2019 (Figure 8A,C), while using measurements from all-three components can lead to identifying 128 sources for the same period (Figure 8B,D).
Given the small number of events, we do not attempt to compute the magnitudes of completeness or -values of the new sources from using different components. However, the exercise shows that using measurements from all-three components can help identify smaller seismic sources (∼M3, Figure 8D).
To compare the method performance of using solely the vertical component [Fan et al., 2018] and combining measurements from all-three components (this study), we use the detected ISC earthquakes that are within the region of interest as benchmark references (Fig-
ure 7). We use their detections to perform a location procedure following the same steps in Section 2.2. The resolved locations and occurrence times are then compared with the ISC origin times and locations (Figure 13). The median distance deviations between the ISC origins and our estimates are 15.8 km, 16.5 km, and 20.3 km for using the vertical, horizontal, and all-three components, respectively (Figure 13). The associated standard deviations (SD, ) are 107.6 km, 118.5 km, and 128.8 km. The median time differences (origin time minus resolved centroid time) are -9.9 s, -6.3 s, and -7.7 s with standard deviations as 99.5 s, 114.2 s, and 110.6 s, respectively. The negative median time differences suggest that the resolved times are later than the origin times, which is expected as we use surface waves to identify the seismic sources instead of using the body waves. The statistics seem to suggest that using the all-three component measures would lead to worse estimates. However, these estimates are obtained from 150,142, and 213 earthquakes for using the vertical, horizontal, and all-three components, respectively.
For a better comparison, we examine the common events that were detected from different measurement sets (Figure 13B,C,E,F). For example, there are 115 common earthquakes detected by using the vertical and horizontal components (Figure 13B,E), and all -9-Confidential manuscript submitted to Geophysical Journal International events detected by using the vertical data are also detected by using the all-three components (150 events, Figure 13C,F). For these 150 events (Figure 13C), the median deviation distance and the associated standard deviation are smaller for using all-three components (14.8 km and 97.3 km) than only using the vertical component (15.8 km and 107.6 km). The resolved event times have similar deviations for using the vertical and and all-three components (Fig-
). The comparisons show that using all-three components can lead to better accuracy and more detections (Figure 7 and 13). Considering the periods of the surface waves used in the study, these uncertainties are within 1/3 wavelengths of the surface waves (∼70 km for 20 s surface waves).
For each resolved seismic source, we require that the event is detected by at least 100 triads to be included in the catalog of robust sources. This requirement is empirically selected and can be adjusted for different applications. The strict requirement assures the reality of the detected signals, but does not address their physical nature. In addition to seismic sources, heterogeneous Earth structure and sharp material interfaces may reflect and convert incoming waves to generate a secondary coherence surface wavefield [Obara and Matsumura, 2010;Maeda et al., 2014;Yu et al., 2017]. Such cases may cause false detections and can be difficult to rule out when only using vertical records and Rayleigh waves. Here by including Love waves, we establish another way to independently examine the sources being seismic events (e.g., Figure 14), as reflection or conversion unlikely produces both coherent Rayleigh and Love wavefields at the same time [e.g., Buehler et al., 2018;Yu et al., 2021].
Our identified sources are valuable as a first-order data product, and we caution that further careful waveform modeling and statistical reasoning are necessary to confirm the physical nature of the events [e.g., Fan et al., 2020Fan et al., , 2022]].
The resolved locations have uncertainties introduced by the array configuration, grid parameterization, and unaccounted arrival angle anomaly that represents the integral of the heterogeneous velocity structures along the ray path [Fan et al., 2018]. One data-driven approach to investigate the location uncertainty is to examine the spatial structure of grids within a pre-defined misfit threshold [de Groot-Hedlin and Hedlin, 2018;Fan et al., 2018].
For example, grids with misfit values within 25% of the minimum value can be considered as possible source locations [Fan et al., 2018], and their distance covariance matrix can be used to obtain an uncertainty ellipse around the location estimate (e.g., red dashline in Figure 14).
The uncertainty ellipse has a minor and major axis, and the major axis length can be used as a metric to describe the the maximum uncertainty of the source location (e.g., the 5 • threshold to determine robust sources). This approach provides one way to formally evaluate the location error for a given cluster of detections, but the misfit threshold is subjectively chosen.
The 25% of the optimal value is a conservative choice, and for a given array configuration, the threshold may be different for sources of different regions.
We can also bootstrap the triads for specific cases to quantify the location uncertainty.
Here we take the 2009 Cascadia Mw 5.7 VLFE as an example to investigate the uncertainties because the case was independently validated by geodetic observations [Fan et al., 2022]. By bootstrapping 1000 times of the subarrays, we can obtain a location uncertainty (standard deviation) of 0.27 • in latitude and 0.26 • in longitude for the example event (Figure 14). Here in each bootstrap run, we randomly draw equal number of triad measurements as used for the location but allow the same measurement being selected multiple times for a realization.
These uncertainty estimates are smaller than the ones obtained from the spatial structure of grids within 25% of the optimal value (maximum uncertainty of 155 km, Figure 14). The bootstrap operation is computationally demanding and might be most suitable for a few selected sources.
The detection method misses some ISC earthquakes because we require a coherent wavefield and the location resolution generally decreases when sources are further away from -10-Confidential manuscript submitted to Geophysical Journal International the seismic arrays [Fan et al., 2018]. As we determine seismic source locations first and then estimate the origin times, the location uncertainty would be translated into timing biases.
Additionally, we simplistically estimate one average surface wave speed for each source. Although the simplification frees us from requiring an Earth velocity model, it also introduce biases. We currently do not account for differences in the Love and Rayleigh wave velocities because mixed surface waves are also used to locate the seismic sources, and it is challenging to asses the type of surface wave dominating the propagation across the triads (e.g., Figure 4, 5, and 14). Future implementation of using arrival times from 3D surface wave tomographic models may help to improve the accuracy of event times.
Our method harnesses the coherence of Rayleigh, Love, and mixed-surface waves. Although having different characteristics, these surface waves share similar ray paths along the great circle paths, enabling us to combine their propagation directions to determine seismic source locations. Such a strategy produces more detections and locations with higher qualities, leading to identifying new seismic sources that may have been missed in standard catalogs [Fan et al., 2019[Fan et al., , 2020[Fan et al., , 2022]]. For example, our method can help detect foreshocks (e.g., Figure 9A), earthquakes in remote regions such as the mid-ocean ridges (e.g., Fig-
ure 9B,C), and offshore earthquakes buried in the noise (e.g., Figure 9D). On July 5, 2019, an earthquake sequence of an 4.7, 5.6, and 4.6 struck near Haida Gwaii in between the Cascadia subduction zone and the Queen Charlotte Fault. We observe an 3.9 ± 0.1 foreshock (magnitude estimate after the scaling correction) from the same region preceding the sequence by about 12 hours (Figure 9A). Similarly, we observe an Mw 4.7 ± 0.1 event near the northern Mid-Atlantic Ridge on July 10, 2017 (Figure 9B) and an Mw 4.4 ± 0.1 near the Reykjanes Ridge on July 22, 2019 (Figure 9C). Intriguingly, we also observe an Mw 3.2 ± 0.2 event near the Camden Bay faults offshore Alaska in the Arctic ocean on July 30, 2019 (Figure 9D). We attribute the detected sources as earthquakes as they are close to well-mapped fault networks. The examples were not reported in the ISC catalog but all generated coherent surface wavefields across continents and can be easily identified from aligned wave records. The cases highlight that our method can aid identifying seismic events in offshore regions. Particularly, seismicity near active spreading centers along mid-ocean ridges and oceanic transform faults are difficult to observe due to the scarcity of offshore instruments [e.g., Parnell-Turner et al., 2022], and our detections in combination with other catalogs may help provide new insights into the processes of crustal deformation near these plate boundaries.
Our algorithm has the advantage in identifying non-earthquake seismic sources [e.g., Fan et al., 2019;Luo et al., 2023]. Particularly, transient surficial processes often produce limited high-frequency body waves but can generate coherent surface waves [e.g., Ekström and Stark, 2013]. Additional, vertical force processes, e.g., volcanic eruptions, may preferably produce either Rayleigh or Love waves but not both types [e.g., Poli and Shapiro, 2022].
In such cases, previous methods may miss these transient cases [Fan et al., 2018]. For example, we located an event with a scaled Mw of 4.7 ± 0.4 in the Gulf of Mexico on July 3, 2019 and an event in Alaska with a Mw of 2.8 ± 0.2 on July 30, 2019. The events are not included in standard catalogs. The seismic source in the Gulf of Mexico may have been a
submarine landslide as suggested in Fan et al. [2020]. The event in Alaska occurred near the Chugach Mountains that are surrounded by seismic instruments. However, no earthquake was reported at the site at our resolved time.
The region had landslides with volumes as great as 15 × 10 6 m 3 [Jibson et al., 2004], and the observed source on July 30, 2019 may have -11-Confidential manuscript submitted to Geophysical Journal International been a landslide with a smaller volume. Further investigation including remote sensing data is warranted for the case.
Glacial quakes were first identified near Greenland [Ekström et al., 2003]. The class of quakes has moderate equivalent earthquake magnitude and occur near major glacier outlets [Nettles and Ekström, 2010]. They are generated by rapid ice motions but with durations one order of magnitude greater than typical tectonic earthquakes [Tsai and Ekström, 2007;Veitch and Nettles, 2012;Olsen and Nettles, 2017]. These events do not produce much high-frequency seismic radiation, but can generate globally recorded coherent surface waves [Nettles and Ekström, 2010]. From 1993 to 2013, there were 444 glacial quakes reported in Greenland [Ekström et al., 2003;Tsai and Ekström, 2007;Veitch and Nettles, 2012;Olsen and Nettles, 2017]. During July 2019, we identified 20 seismic sources in Greenland with the majority of them near Jakobshavn Isbrae, Daugaard Jensen, and northwest Greenland (Figure 11A). These events generated globally recorded coherent surface wavefields (e.g., Figure 11B). We speculate that these events are likely glacial quakes. The events all have scaled magnitudes ( ) closely around 4.5 (Figure 11C), which is similar to the characteristic magnitude reported in previous studies [Nettles and Ekström, 2010]. We found that the events occurred intermittently during the study period and that multiple cases may occur as clusters with short inter-event intervals (Figure 11C). These examples suggest that combining
Love and Rayleigh waves is promising in identifying more exotic seismic sources such as the glacial quakes in Greenland.
We have developed an approach to combine Love and Rayleigh waves from records of all three components for detecting and locating seismic sources. The new method uses intermediate-period (20-50 s) surface waves and globally distributed arrays that do not need to be spatially adjacent. By taking advantage of local coherence, the method measures surface wave propagation directions and centroid arrival times, and these measurements are aggregated to locate seismic sources. The method assumes that surface waves propagate along great circle paths with small deviations, and that Rayleigh and Love waves share similar ray paths. We do not require much a priori information on the seismic events or a velocity model to locate the sources. The method is one of the few to incorporate Love waves in identifying seismic sources, and the combination of both Love and Rayleigh waves can enhance detection sensitivity, superior than only using Rayleigh waves from the vertical component.
Taking July 2019 as an example month, we show that the method is capable of detecting earthquakes as well as exotic sources. The method is particularly useful in identifying remote events at mid-ocean ridges or Greenland that are missed in standard catalogs. Our results also suggest a new direction for examining data from other geophysical instruments, such as the tiltmeters of the Japanese Hi-net array [Okada et al., 2004].
Management Center (DMC) of the Incorporated Research Institutions for Seismology (IRIS)
and European Integrated Data Archive (EIDA). Seismic networks used in this study are listed in an electronic supplement with their DOIs included. The seismic data were downloaded using ObsPy [Beyreuther et al., 2010], and the synthetic seismograms were computed using Instaseis [Driel et al., 2015]. The AELUMA code can be obtained on request through IRIS [Hutko et al., 2017]. The seismic source catalog of July 2019 from this study is included as an electronic supplement. -14-Confidential manuscript submitted to Geophysical Journal International Figure 2. Globally distributed searching points (A) and example dense points in California (B). Candidate locations are denoted as red dots. Plate boundaries are marked as blue lines in (A), and local faults are marked
, and
) independently or collectively.2.3 Seismic Source MagnitudeOnce a source location is determined, the source magnitude M SE can be estimated using the measured beampower values. For each measurement that is used for the final location procedure, an estimate can be computed as M SE = log 10 + 1.66 log 10 + 2, (2)