https://orcid.org/0000-0002-2983-8240
SUMMARY
An earthquake rupture process can be kinematically described by rupture velocity, duration and spatial extent. These key kinematic source parameters provide important constraints on earthquake physics and rupture dynamics. In particular, core questions in earthquake science can be addressed once these properties of small earthquakes are well resolved. However, these parameters of small earthquakes are poorly understood, often limited by available data sets and methodologies. The Incorporated Research Institutions for Seismology Community Wavefield Experiment in Oklahoma deployed ∼350 three-component nodal stations within 40 km2 for a month, offering an unprecedented opportunity to test new methodologies for resolving small earthquake finite source properties in high resolution. In this study, we demonstrate the power of the nodal data set to resolve the variations in the seismic wavefield over the focal sphere due to the finite source attributes of an M2 earthquake within the array. The dense coverage allows us to tightly constrain rupture area using the second moment method even for such a small earthquake. The M2 earthquake was a strike-slip event and unilaterally propagated towards the surface at 90 per cent local S-wave speed (2.93 km s−1). The earthquake lasted ∼0.019 s and ruptured Lc ∼70 m and Wc ∼45 m. With the resolved rupture area, the stress-drop of the earthquake is estimated as 7.3 MPa for Mw 2.3. We demonstrate that the maximum and minimum bounds on rupture area are within a factor of two, much lower than typical stress-drop uncertainty, despite a suboptimal station distribution. The rupture properties suggest that there is little difference between the M2 Oklahoma earthquake and typical large earthquakes. The new three-component nodal systems have great potential for improving the resolution of studies of earthquake source properties.
1 INTRODUCTION
Rupture velocity, duration and spatial extents directly characterize basic earthquake rupture processes. These key kinematic finite source parameters can provide observational constraints for earthquake dynamic processes. For instance, earthquake stress-drop, which is closely related to near-field ground motion, can be directly estimated from the earthquake moment and rupture area (Brune 1970; Madariaga 1976). Therefore, accurate, model-free estimates of these key kinematic parameters will help in understanding earthquake physics, fault zone properties and the conditions of rupture dynamics (e.g. Beroza & Ellsworth 1996; Mai & Beroza 2000; Archuleta & Ji 2016; Meier et al.2017; Melgar & Hayes 2017; Thingbaijam et al.2017). However, even for moderate earthquakes, these parameters sometimes are poorly constrained because of the limitations of available data sets and current analysis methods (e.g. Venkataraman et al.2000). For small earthquakes, these kinematic attributes are often assumed to be unresolveable and the events are represented with simple models despite near-field observations suggesting the opposite (e.g. Lanza et al.1999; Ide 2001; McGuire 2004; Yamada et al.2005; McGarr et al.2009; Boettcher et al.2015). Moderate and small earthquake rupture processes are often examined through spectral fitting techniques, which require an a priori physical rupture model (e.g. Abercrombie 1995; Imanishi et al.2004; Shearer et al.2006; Noda et al.2013). However, numerical experiments by Kaneko & Shearer (2014, 2015) demonstrated that standard spectral fitting can lead to roughly one order of magnitude variation in stress-drop estimates that do not reflect the actual rupture properties even for the simple crack model. These data and methodology limitations make some of the fundamental earthquake questions elusive. For example, do small and large earthquakes rupture in similar fashions (e.g. Dreger et al.2007; Meier et al.2016)? Is the large scatter in stress-drop estimates due to earthquake complexities or caveats of the spectral fitting techniques (e.g. Allmann & Shearer 2009, 2007)? To address these questions, it is important to have dense near-source seismograph instrumentation and kinematic finite source measurements that are independent from assumed rupture models (McGuire 2004, 2017).
Here, we apply the second moment method to resolve the rupture process of one local M2 earthquake recorded by the Incorporated Research Institutions for Seismology (IRIS) Community Wavefield Demonstration Experiment in Oklahoma (Anderson et al.2016; Sweet et al.2018). The second-degree moments of an earthquake describe the length, width, duration and directivity of its rupture (Backus & Mulcahy 1976a,b; Backus 1977a,b; Doornbos 1982a,b; Stump & Johnson 1982; Silver 1983; Gusev & Pavlov 1988; Bukchin 1995; Das & Kostrov 1997; McGuire et al.2001; Clévédé et al.2004). The second moments can be inverted for using measurements of apparent duration of any particular seismic phase at different stations for a given earthquake. The apparent durations are measured from the apparent source time functions (ASTF) for each station and seismic phase by removing the path effects, which is often done by empirical Green’s function (EGF) deconvolution (McGuire 2004, 2017). The second moment approach does not require an a priori physical rupture model, and is particularly useful to resolve moderate to small earthquake (Mw ≤ 6) rupture processes when near-field geodetic observations are absent (e.g. McGuire 2004; Chen & McGuire 2016; Gong & McGuire 2018).
Aiming to observe the seismic wavefield with minimal aliasing and maximum resolution, IRIS conducted a wavefield experiment at north central Oklahoma (OK) in summer 2016 (Anderson et al.2016; Sweet et al.2018). The experiment installed 386 seismographs, including 359 three-component short-period seismometres (nodes), 18 broad-band stations and 9 infrasound stations (Fig. 1). The stations were deployed within ∼5 km × ∼8 km with the nodal array recording about one month of data. The nodal array captured at least two M2 earthquakes during its deployment (Fig. 1). Compared to typical data sets (∼20 observations), the wavefield experiment offers an unprecedented opportunity to evaluate small earthquake finite source attributes and a platform to test new finite source estimation techniques that take advantage of the next-generation seismographs.
Station distribution of the IRIS Community Wavefield Demonstration Experiment in Oklahoma. The coloured lines show three seismic line profiles and the grey squares are the rest of the array. The 7/11 M2 earthquake of interests is labelled as a red cross with its focal mechanism by the side. The dense instrumentation yields high location precision with horizontal uncertainties within 100 m. The upper left inset shows the array location in central northern Oklahoma (red square). The lower right inset shows the local 1-D velocity model with the basement of the sedimentary layer set to be 2270 m.
2 DATA: IRIS COMMUNITY WAVEFIELD DEMONSTRATION EXPERIMENT IN OKLAHOMA
The Oklahoma wavefield experiment had four main components, including 3 seismic lines, a 7-layer gradiometer, 18 broad-band stations and 9 infrasound stations (Fig. 1). The seismic lines and nested gradiometer are nodal sensors. The experiment array geometry was determined by permitting conditions in the field. In this study, we focus on the three seismic lines. The three seismic lines are composed of 220 5-Hz three-component nodal stations with interstation spacing ∼100 m following roads. These three lines were initially designed in anticipation of a controlled source experiment. During the wavefield experiment, the nodal stations record signals at a 250 Hz sampling rate. The experiment field site has a relatively simple velocity structure (Schoenball & Ellsworth 2017). Our following analysis is based on the 1-D velocity model used in Schoenball & Ellsworth (2017) and by the Oklahoma Geological Survey (OGS; Darold et al.2015) with the top of basement set as 2270 m (adjustment suggested by Maria Beatrice Magnani, personal communication, 2017; inset of Fig. 1).
Two M2 earthquakes occurred beneath the array during the nodal array deployment (Fig. 1). The 2016 July 15 M2 earthquake had only two usable EGFs (consecutively ruptured within 5 s) for second moments analysis. Here, the useable EGFs can be clearly visually identified with sharp onsets, have high signal-to-noise ratios (SNRs; root mean squares of 0.3 s signal windows are at least three times larger than those of the 0.3 s preceding time window) and share similar moveout pattern as the main shock. To identify the moveout pattern, we apply cross-correlation to empirically align the main-shock records (Houser et al.2008) and then apply the obtained alignment to continuous 1-hr records after the main-shock origin time. EGFs close to the main shock would have a vertical moveout pattern, while events away from the main shock would have a skewed moveout pattern. Due to uneven temporal distribution of these EGFs, there is no temporal preference. To minimize the EGF finite source effects, EGFs with small amplitudes but low noise levels are preferred.
We focus on the source process of the 2016 July 11 M2 earthquake, which was hosted by an active fault with strike at ∼240°. We first bandpass three components of the stations at 10–100 Hz with a second-order Butterworth filter, then cross-correlate the records to pick P- and S-wave arrivals (Houser et al.2008). With our preferred velocity model, 3-D grid search is performed to fit the S–Ptimes at each station of the three line profiles. The M2 earthquake is located at latitude 36.617°, longitude −97.687° at 2.8 km depth (Table 1). We then visually identify the P-wave polarities at each station to resolve the focal mechanism of the earthquake with HASH (Hardebeck & Shearer 2002, 2003). The preferred focal mechanisms are 328.0°/80.0°/−176.0° and 237.3°/86.1°/−10° for strike, dip and rake. The 237.3°/86.1°/−10° fault plane is preferred because of its close agreement with the fault strike delineated by local seismicity (catalogue provided by Heather DeShon and Louis Quinones, personal communication, 2018; Fig. 1).
M2 earthquake and its second moments.
| Time | Location | Fault geometry | Magnitude | |
|---|---|---|---|---|
| (UTC) | [Lat°/lon°/depth (km)] | (Strike/dip/rake) | (ML) | |
| 2016/7/11 5:55:16.5 | 36.617/−97.687/2.8 | 237.3/86.1/−10 | 2.3 | |
| Length (Lc, m) | Width (Wc, m) | Vs (km s−1) | Vd (km s−1) | Duration (s) |
| 71.2 | 44.5 | −0.58 | −2.87 | 0.0192 |
| Time | Location | Fault geometry | Magnitude | |
|---|---|---|---|---|
| (UTC) | [Lat°/lon°/depth (km)] | (Strike/dip/rake) | (ML) | |
| 2016/7/11 5:55:16.5 | 36.617/−97.687/2.8 | 237.3/86.1/−10 | 2.3 | |
| Length (Lc, m) | Width (Wc, m) | Vs (km s−1) | Vd (km s−1) | Duration (s) |
| 71.2 | 44.5 | −0.58 | −2.87 | 0.0192 |
Vs: rupture velocity along strike, Vd: rupture velocity along dip and ML: local magnitude.
M2 earthquake and its second moments.
| Time | Location | Fault geometry | Magnitude | |
|---|---|---|---|---|
| (UTC) | [Lat°/lon°/depth (km)] | (Strike/dip/rake) | (ML) | |
| 2016/7/11 5:55:16.5 | 36.617/−97.687/2.8 | 237.3/86.1/−10 | 2.3 | |
| Length (Lc, m) | Width (Wc, m) | Vs (km s−1) | Vd (km s−1) | Duration (s) |
| 71.2 | 44.5 | −0.58 | −2.87 | 0.0192 |
| Time | Location | Fault geometry | Magnitude | |
|---|---|---|---|---|
| (UTC) | [Lat°/lon°/depth (km)] | (Strike/dip/rake) | (ML) | |
| 2016/7/11 5:55:16.5 | 36.617/−97.687/2.8 | 237.3/86.1/−10 | 2.3 | |
| Length (Lc, m) | Width (Wc, m) | Vs (km s−1) | Vd (km s−1) | Duration (s) |
| 71.2 | 44.5 | −0.58 | −2.87 | 0.0192 |
Vs: rupture velocity along strike, Vd: rupture velocity along dip and ML: local magnitude.
With the earthquake location, three components at each station are rotated into vertical, radial and transverse components (Fig. S1–S3, Supporting Information). As shown in Fig. S1 of the Supporting Information, the P waves are complex with multiple arrivals in a short time window. Therefore, we do not use P waves for second moment analysis because of the interference between the direct P waves and these high amplitude phases (e.g. multiples and converted phases). We focus on using the transverse component (SH) S waves to resolve the ASTFs and do not use the radial component because of redundant information and lower SNRs. To solve the ASTFs at each station, we use seven EGFs that are detected within 1 hr after the M2 main shock. These EGFs are recorded with peak S-wave amplitudes at least 100 times smaller than the main shock at a given station (Fig. 2; Figs S4 and S5, Supporting Information). SH waves of the main shock and the EGFs are tapered into 0.3 s windows by a Tukey window with the first and last 2.5 per cent samples equal to parts of a cosine. At each station, a given EGF is cross-correlated with the main shock and is shifted 0.04 s ahead to assure causality for the later resolved ASTFs.
Example of SH apparent duration measurements. Each column (A–D) shows results of a different station along Line 1. The first row shows the data fit (records are self-normalized to unity). The second row compares one EGF with the main-shock records (records are self-normalized to unity). The third row shows the misfit curves of both the ending point and the starting point. The fourth row shows the obtained apparent source time functions with the red dots indicating centroid times and the red bars indicating the apparent durations.
3 METHODS
3.1 Second moments inversion
Example of SH data fit and resolved apparent source time functions (ASTF) with one EGF for three seismic line profiles.
Apparent duration (|$\tau _c(\underline{s})$|) measurements from seven EGFs of the three seismic line profiles. The first row shows the measurements obtained from EGF deconvolution. Different colours represent different EGFs. The second row shows all the measurements and cubic smoothing splines fitted to the measurements. The third row shows the data after removing the outliers and their fitted spline functions of each line profiles.
For a given station, the estimates of apparent duration, |$\tau _c(\underline{s})=2\sqrt{ \mu ^{(0,2)}(\underline{{s}})}$|, from different EGFs generally agree with each other, albeit outliers with extreme values may be present (Fig. 4). We first conservatively remove measurements with |$\tau _c(\underline{s})$| less than 0.008 s (two sample points). As we know |$\tau _c (\underline{s})$| varies smoothly over space (McGuire 2004, 2017), this useful a priori information can be used to efficiently remove the outliers. For instance, we first combine all the measurements for Line 1 and then fit a cubic smoothing spline to the data set (Fig. 4B). Each data point is inversely weighted by its misfit and a smoothing parameter is applied to obtain the cubic smoothing spline. The smoothing parameter is calculated from |$p=\frac{1}{1+h^3/6}$|, where h is the interstation spacing (0.1 km; Thompson & Shure 2016). When p is set to 0, a straight line is fitted to the data set and p > 1 will lead to a rougher variational cubic spline interpolant. After obtaining the smooth spline, we remove data points that are one standard deviation (SD) away from the spline, where the SD is calculated from all the measurements (Fig. 4C). A new cubic smoothing spline function is fitted for the remaining data set for visualization (Fig. 5). All the remaining data (658 measurements) are used to invert the second moments for the M2 earthquake.
Map view of the cubic smoothing splines from the cleaned data (Fig. 4C). The lower left inset shows the measurement distribution. The lower right inset is the polar plot of cubic smoothing splines, in which radii are the incident angles at each stations and azimuths are calculated with respect to the M2 main-shock epicentre.
3.2 Uncertainty: estimating upper and lower bounds on rupture area
3.3 Corner frequency (fc) and fall-off exponent (n)
Earthquake source spectra are often modelled to investigate the finite source attributes of earthquakes (e.g. Houston & Kanamori 1986; Houston 1990; Prieto et al.2004; Baltay et al.2010; Chen & Shearer 2011, 2013; Baltay et al.2013; Uchide et al.2014; Goebel et al.2015; Imanishi & Uchide 2017). Two key parameters describing an earthquake spectrum are the average corner frequency (fc) and fall-off exponent (n). The average corner frequency can be simply linked to earthquake static stress-drop (Δσ) once a dynamic rupture model is assumed (Eshelby 1957; Brune 1970; Sato & Hirasawa 1973; Madariaga 1976; Boatwright 1980).
Here we use the spectral ratio method to estimate S-wave corner frequencies and fall-off exponents at each stations (e.g. Ide et al.2003; Imanishi et al.2004; Abercrombie 2014; Lengliné et al.2014; Uchide et al.2014; Huang et al.2016, 2017). We use the same seven EGFs used for the second moment inversion to obtain the spectral ratios of the M2 earthquake at each station. With the spectral ratios, we then solve for the corner frequencies and fall-off exponents.
In practice, at a given station, we first remove the path effects by deconvolving an EGF S-wave from the M2 main-shock record in the frequency domain to obtain a source spectrum. The S-wave window is 2 s long, which is long enough for obtaining stable solutions. A 5 per cent Tukey window is applied before calculating the spectrum. The EGF deconvolution is performed with all the seven EGFs at both transverse and radial components separately. We then derive the S-wave source spectrum at a given station by geometrically averaging all the obtained source spectra of different stations within 1 km of that station. The geometric mean is performed here to have the arithmetic mean of the spectra in the log space. Finally, we apply a nonlinear least-squares algorithm to fit the spectra from ∼5 to ∼35 Hz with the Boatwright spectral model assuming γ = 2 (Boatwright 1980; Huang et al.2017; Imanishi & Uchide 2017). The spectra fitting range is suggested in Huang et al. (2017) for M2-induced earthquakes in Arkansas. No other smoothing or post-processing is applied to the spectra.
4 RESULTS
4.1 Second moments
The estimates of the second moments are listed in Table 1. The results suggest that the M2 earthquake ruptured unilaterally with optimal estimates of Lc as 71.2 m and Wc as 44.5 m. The rupture propagated near vertically towards the surface at ∼90 per cent local shear wave speed (|v0|=2.93 km s−1) and had a characteristic duration of τc = 0.019 s. The horizontal component of rupture propagation was to the northeast (−0.58 km s−1). The resolved M2 earthquake finite source model can explain the observations well with residuals below the sampling rate (Fig. 6). As shown in Fig. 6, predictions of |$\tau _c(\underline{s})$| at each station agree well with the observations and share resemblances with the interpolated cubic smoothing splines of the three seismic lines despite no requirement that they have similar functional forms. With the inverted rupture length and width, assuming the local magnitude is the same as the moment magnitude (Sumy et al.2017), Mw2.3 (M0 = 3.16 × 1012 Nm), the stress-drop of the earthquake is estimated as 7.3 MPa using the solutions for an elliptical crack. Due to limited observations, a scaling relationship between ML and Mw is not well defined for the wavefield experiment region. To account for the magnitude uncertainties, we also calculate stress-drops with Mw 2.2 and Mw 2.4 for the M2 main shock (Table 2).
Predictions of apparent durations (|$\tau _c(\underline{s})$|) at each station from the inverted finite source model (Table 1). (A) Map view of the predictions, the legends are similar with Fig. 5 except the lower left inset showing the residual distribution. (B–D) Data fits to the apparent duration measurements of three seismic line profiles.
Stress-drop estimates (Δσ MPa) from second moments.
| Length (Lc, m) | Width (Wc, m) | Mw = 2.3 | Mw = 2.2 | Mw = 2.4 | |
|---|---|---|---|---|---|
| Δσ0 | 71.2 | 44.5 | 7.3 | 5.2 | 10.34 |
| |$\Delta \sigma _\mathrm{ d}^{\mathrm{ L}}$| | 78.3 | 46.3 | 6.2 | 4.4 | 8.7 |
| |$\Delta \sigma _\mathrm{ d}^{\mathrm{ U}}$| | 63.5 | 40.5 | 9.9 | 7.0 | 14.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ sL}}$| | 93.1 | 47.0 | 5.0 | 3.5 | 7.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ sU}}$| | 74.1 | 39.1 | 9.1 | 6.4 | 12.8 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ aL}}$| | 94.4 | 50.6 | 4.3 | 3.0 | 6.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ aU}}$| | 43.1 | 21.0 | 54.0 | 38.2 | 76.2 |
| Length (Lc, m) | Width (Wc, m) | Mw = 2.3 | Mw = 2.2 | Mw = 2.4 | |
|---|---|---|---|---|---|
| Δσ0 | 71.2 | 44.5 | 7.3 | 5.2 | 10.34 |
| |$\Delta \sigma _\mathrm{ d}^{\mathrm{ L}}$| | 78.3 | 46.3 | 6.2 | 4.4 | 8.7 |
| |$\Delta \sigma _\mathrm{ d}^{\mathrm{ U}}$| | 63.5 | 40.5 | 9.9 | 7.0 | 14.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ sL}}$| | 93.1 | 47.0 | 5.0 | 3.5 | 7.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ sU}}$| | 74.1 | 39.1 | 9.1 | 6.4 | 12.8 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ aL}}$| | 94.4 | 50.6 | 4.3 | 3.0 | 6.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ aU}}$| | 43.1 | 21.0 | 54.0 | 38.2 | 76.2 |
Δσ0 is the optimal stress-drop estimate, |$\Delta \sigma _\mathrm{ d}^{\mathrm{ L}}$| and |$\Delta \sigma _\mathrm{ d}^{\mathrm{ U}}$| are the stress-drop lower and upper bounds derived from measurement bootstrapping, |$\Delta \sigma _{\mathrm{ m}}^{\mathrm{ sL}}$| and |$\Delta \sigma _{\mathrm{ m}}^{\mathrm{ sU}}$| are the stress-drop lower and upper bounds derived from maximizing rupture area based on 95 per cent χ2 confidence level by fitting the cubic smoothing splines, |$\Delta \sigma _{\mathrm{ m}}^{\mathrm{ aL}}$| and |$\Delta \sigma _{\mathrm{ m}}^{\mathrm{ aU}}$| are the stress-drop lower and upper bounds derived from maximizing rupture area based on 95 per cent χ2 confidence level by fitting all the data.
Stress-drop estimates (Δσ MPa) from second moments.
| Length (Lc, m) | Width (Wc, m) | Mw = 2.3 | Mw = 2.2 | Mw = 2.4 | |
|---|---|---|---|---|---|
| Δσ0 | 71.2 | 44.5 | 7.3 | 5.2 | 10.34 |
| |$\Delta \sigma _\mathrm{ d}^{\mathrm{ L}}$| | 78.3 | 46.3 | 6.2 | 4.4 | 8.7 |
| |$\Delta \sigma _\mathrm{ d}^{\mathrm{ U}}$| | 63.5 | 40.5 | 9.9 | 7.0 | 14.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ sL}}$| | 93.1 | 47.0 | 5.0 | 3.5 | 7.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ sU}}$| | 74.1 | 39.1 | 9.1 | 6.4 | 12.8 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ aL}}$| | 94.4 | 50.6 | 4.3 | 3.0 | 6.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ aU}}$| | 43.1 | 21.0 | 54.0 | 38.2 | 76.2 |
| Length (Lc, m) | Width (Wc, m) | Mw = 2.3 | Mw = 2.2 | Mw = 2.4 | |
|---|---|---|---|---|---|
| Δσ0 | 71.2 | 44.5 | 7.3 | 5.2 | 10.34 |
| |$\Delta \sigma _\mathrm{ d}^{\mathrm{ L}}$| | 78.3 | 46.3 | 6.2 | 4.4 | 8.7 |
| |$\Delta \sigma _\mathrm{ d}^{\mathrm{ U}}$| | 63.5 | 40.5 | 9.9 | 7.0 | 14.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ sL}}$| | 93.1 | 47.0 | 5.0 | 3.5 | 7.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ sU}}$| | 74.1 | 39.1 | 9.1 | 6.4 | 12.8 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ aL}}$| | 94.4 | 50.6 | 4.3 | 3.0 | 6.0 |
| |$\Delta \sigma _\mathrm{ m}^{\mathrm{ aU}}$| | 43.1 | 21.0 | 54.0 | 38.2 | 76.2 |
Δσ0 is the optimal stress-drop estimate, |$\Delta \sigma _\mathrm{ d}^{\mathrm{ L}}$| and |$\Delta \sigma _\mathrm{ d}^{\mathrm{ U}}$| are the stress-drop lower and upper bounds derived from measurement bootstrapping, |$\Delta \sigma _{\mathrm{ m}}^{\mathrm{ sL}}$| and |$\Delta \sigma _{\mathrm{ m}}^{\mathrm{ sU}}$| are the stress-drop lower and upper bounds derived from maximizing rupture area based on 95 per cent χ2 confidence level by fitting the cubic smoothing splines, |$\Delta \sigma _{\mathrm{ m}}^{\mathrm{ aL}}$| and |$\Delta \sigma _{\mathrm{ m}}^{\mathrm{ aU}}$| are the stress-drop lower and upper bounds derived from maximizing rupture area based on 95 per cent χ2 confidence level by fitting all the data.
We first evaluate the uncertainties of the resolved second moments with bootstrap resampling of the measurements (Efron & Tibshirani 1994). We randomly draw 50 per cent of all the measurements (329 out of 658 measurements), which may not sample all three seismic lines, and then invert for the second moments of the earthquake following the same method described above (Fig. 7). The procedure is performed 5000 times for obtaining distributions of the second moments. The results show that the resolved second moments are tightly bounded with small SDs, suggesting that the earthquake finite source attributes are well resolved. The uncertainties are further evaluated with bootstrapping all the measurements including outliers that were removed from the cubic smoothing spline procedure. The results show that the means of the second moments are well constrained, while the SDs would increase when extreme outliers are present (Fig. S6, Supporting Information). Considering within one SD of the second moments in Fig. 7, the stress-drop of the earthquake varies from 6.2 to 9.9 MPa (Table 2).
Distributions of characteristic length (A), width (B), duration (C) and rupture velocity (D) obtained by bootstrap resampling. μ and σ are the mean and SD for the parameters.
Inaccurate earthquake location or focal mechanism can introduce uncertainties in the resolved second moments. Earthquake location and focal mechanism have direct influences on |$\underline{\underline{A}}$| in eq. (8), but do not impact the measurements (|$\underline{b}$| in eq. 8). Because of the dense instrumentation, the earthquake horizontal location can be well resolved within 100 m uncertainty (Figs S1 and S2, Supporting Information). The resolution of absolute earthquake depth is largely determined by the velocity model (e.g. Lin et al.2007). In Oklahoma, the absolute vertical precision peaks around 500 m and vertical precisions of 77 per cent of events are within 1 km (Schoenball & Ellsworth 2017). To understand the uncertainties introduced by inaccurate earthquake location, we perform second moments inversion with source location grids within [−100 100] m along latitude by [−100 100] m along longitude by [−500 500] m in depth of the preferred location (36.617°/ − 97.687°/2.8 km). The grids are separated 50 m in distance, leading to 525 grid points in total. As shown in Figs 8(A)–(E), the means and SDs of the second moments are similar to those derived from bootstrap analysis (Fig. 7). Because of the good azimuthal station coverage and clear onsets of the P-wave first motions, the focal mechanism is uniquely determined by HASH for this case. To understand the influences of possible inaccurate focal mechanism (strike and dip), we invert the second moments with strike and dip varying from −10° to −10° of the preferred fault geometry (strike: 237.3°, dip: 86.1°) with 1° searching step. Similarly, the resolved means and SDs of the second moments do not vary much from those obtained from bootstrap analysis (Figs 8F–J). Consequentially, the stress-drop estimates are well-resolved given the range of varying source locations and focal mechanisms (6.0–10.6 MPa). The results show that the obtained finite source attributes of the earthquake are robust.
Distributions of characteristic length, width, duration, rupture velocity and stress-drop caused by uncertainties in earthquake location (first row) and focal mechanism (second row). μ and σ are the mean and SD for the parameters.
Following the method described in Section 3.2, we finally estimate maximum/minimum rupture area that is permitted by the data set within 95 per cent confidence level. We then use the estimated upper and lower rupture area to constrain the minimum and maximum possible stress-drops. These stress-drop limits are different from the ones derived from the bootstrap analysis above and are perhaps more accurate as they specifically search for the extreme values given the full data set. For robust estimations, we use interpolated cubic smoothing splines of three seismic lines as input to stabilize the inversion. The stress-drop of the earthquake (Mw 2.3) varies from 5.0 MPa for maximum rupture area to 9.1 MPa for minimum rupture area (Fig. 9 and Table 2). Due to uneven data distributions of three seismic profiles, the estimated stress-drop varies from 4.3 MPa for maximum rupture area to 54.0 MPa for minimum area when all the data are used for inversion. It is worth to note that rupture velocity and rupture area are directly related because of the well-resolved rupture duration (τc ≃ 0.02 s): maximizing the rupture area will increase the rupture velocity and minimizing the rupture area will slow down the rupture propagation. The rupture velocities corresponding to maximum/minimum rupture area range from 3.8 to 2.88 km s−1 when fitting the cubic smoothing splines. The resolved minimum area finite source model with all the data suggests a very slow rupture, 1.27 km s−1, which leads to the rupture length and width less than half of these parameters of the preferred model. As shown in Fig. 9(C), the largest prediction discrepancies between this extreme model to the rest would be present at the western edge of Line 1. Unfortunately, we do not have high SNR records in that region for a verdict. The predictions of the models fall within the measurements spread, albeit less well fitted at the northern end of Line 2 (Fig. 9). The resolved maximum and minimum rupture areas, and hence stress-drops, represent two extreme rupture scenarios that are permitted by the data (Fig. 9 and Table 2).
Predictions of apparent durations at each station from the inverted finite source model with maximum or minimum rupture area (Table 2). (A) Map view of the predictions of inverted source model with maximum rupture area by fitting the cubic smoothing spline, the legends are similar to Fig. 6. (B) Map view of the predictions of inverted source model with minimum rupture area by fitting the cubic smoothing spline. (C–E) Data fits to the apparent duration measurements of three seismic line profiles. (C)–(E) share the same legend as in (C). Predictions with inverted models with maximum/minimum rupture area by fitting all the data are also included for comparison.
4.2 Source spectra
The source spectra (corner frequency and fall-off exponent) show strong directivity effects (Figs 10 and 11). Stations in the rupture direction have corner frequencies as high as 23.2 Hz, while stations away from rupture direction have lower corner frequencies (Fig. 10). The corner frequencies vary with incident angles corresponding to the rupture direction, which inversely correlate to the apparent duration predicted from the resolved finite source model at each station. On the other hand, the fall-off exponent has more complex spatial pattern than the corner frequency (Fig. 11). Source spectra away from the rupture direction show ‘double humps’, suggesting multiple subevents comprise the M2 main shock (e.g. Stations A, B and D in Fig. 10). This observation indicates that the M2 earthquake had a more complicated rupture evolution than a theoretical crack rupture model. Due to different EGFs and different available records for averaging (Figs S4 and S5, Supporting Information), the relative moments are slightly different for each station as shown in Fig. 10. Normalizing the relative moments does not change the resolved corner frequency and fall-off exponent. Without knowing the EGF moments, for simplicity, we keep the relative moments as shown in Fig. 10. The overall averaged spectrum is smooth and can be well fitted with the Boatwright model with corner frequency of 17.1 Hz and fall-off exponent of 2. The averaged spectrum is different from individual spectra and the SD of the measured corner frequencies is 2.3 Hz among all the stations with maximum corner frequency fmax = 23.8 Hz and minimum corner frequency fmin = 12.0 Hz (Fig. 11).
Example S-wave spectra of three stations and the overall average spectra. The relative moments might be different due to different stations having different EGFs and the number of available records for averaging. Spectra of Station A–D are the geometric mean derived from 133, 200, 191 and 203 spectra separately. The overall average spectrum is derived from 2164 spectra.
Map view of the measured corner frequencies (A) and fall-off exponents (B). Black squares show the example stations in Fig. 10.
5 DISCUSSIONS
Earthquake stress-drop provides direct constraints on earthquake scaling and insights of the hosting fault environments (e.g. Abercrombie 1995; Abercrombie & Rice 2005; Chen & Shearer 2013; Cotton et al.2013; Bilek et al.2016; Thingbaijam et al.2017). Therefore, stress-drop from source spectral fitting techniques has been extensively studied at both regional and global scales (e.g. Allmann & Shearer 2007; Denolle & Shearer 2016). Previously, anthropogenic-induced earthquakes (e.g. Keranen et al.2013) are reported to have possible lower stress-drops compared to tectonic earthquakes in some studies (e.g. Goertz-Allmann et al.2011; Hough 2014; Sumy et al.2017; Trugman et al.2017), while the others suggest no significant differences between the two types (e.g. Huang et al.2016, 2017), highlighting the importance of accurate estimates of stress-drops. However, accurate determination of stress-drops is challenging because of uncertainties from both the data and incomplete knowledge of the dynamic rupture models (Shearer et al.2006; Kaneko & Shearer 2014, 2015). In particular, Kaneko & Shearer (2014, 2015) pointed out that corner frequency is largely dependent on take-off angle relative to the source, which can result in a very complex spatial distribution and cause at least one order of magnitude uncertainty in stress-drop estimates. The dense instrumentation of the Oklahoma wavefield experiment provides a great opportunity to investigate the complex spatial pattern of corner frequency and the associated uncertainties on stress-drop estimates.
The stress-drop from the overall averaged corner frequency is 11.4 MPa for κ = 0.26 and 9.1 MPa for κ = 0.28 (Table 3). The κ values are suggested for asymmetric elliptical (κ = 0.26) and circular (κ = 0.28) rupture scenarios in Kaneko & Shearer (2014, 2015). The estimates are slightly higher than the stress-drop estimated from the second moments (7.3 MPa). This is expected as the stations are all in the rupture direction hemisphere, which would lead to a higher averaged corner frequency than the whole sphere average. The general agreement between stress-drops of the corner frequency and the second moment methods is probably because the M2 earthquake ruptured at 90 per cent local shearer wave speed, which is the assumed rupture velocity for corner-frequency-estimated stress-drops. However uncertainties in κ (Kaneko & Shearer 2014, 2015) may lead to one magnitude variation in the stress-drop estimates (Table 3). In addition, due to the complex spatial pattern of corner frequencies (Fig. 11) and large span of measured corner frequencies (12.0 to 23.8 Hz), stress-drop estimated from a small number of seismic stations may be less well constrained as limited spectra measurements might introduce certain biases comparing to the estimate of spherically averaged corner frequency (Kaneko & Shearer 2014, 2015). Comparatively, the stress-drop estimated from the second moments is well bounded from ∼5 to ∼10 MPa. The lower limit (∼5 MPa) is tightly resolved from both data resampling and rupture area uncertainty analyses (both are data driven without a priori assumptions, Table 2), and the upper limit is well resolved around ∼10 MPa, except when using all the data to estimate the minimum area. For the rupture area variation related stress-drop uncertainties, we prefer the 5.0 to 9.1 MPa bound. The preference is chosen not because of its smaller range, but because the cubic smooth splines sample observations at all azimuth evenly comparing to using all the measurements (4.3 to 54.0 MPa). For example, the all-data-minimum-area model predicts significantly smaller apparent durations (|$\tau _c(\underline{s})$|) at the western end of Line 1, which is different from the predictions of the rest models (Fig. 9C).
Stress-drop estimates (Δσ MPa) from corner frequency.
| Mw = 2.3 | Mw = 2.2 | Mw = 2.4 | |
|---|---|---|---|
| κ = 0.21 (Brune 1970) | |$21.6^{+10.1 }_{ -7.7}$| | |$15.3^{+7.2 }_{ -5.5}$| | |$30.5^{+14.3 }_{ -10.9}$| |
| κ = 0.26 (Kaneko & Shearer 2014, 2015) | |$11.4^{+5.3 }_{ -4.1}$| | |$8.0^{+3.8 }_{ -2.9}$| | |$16.0^{+7.5 }_{ -5.7}$| |
| κ = 0.28 (Kaneko & Shearer 2014, 2015) | |$9.1^{+4.27 }_{ -3.25}$| | |$6.4^{+3.0 }_{ -2.3}$| | |$12.8^{+6.0 }_{ -4.6}$| |
| κ = 0.32 (Huang et al.2017) | |$6.1^{+2.9 }_{ -2.2}$| | |$4.3^{+2.0 }_{ -1.5}$| | |$8.6^{+4.0 }_{ -3.1}$| |
| κ = 0.372 (Madariaga 1976) | |$3.9^{+1.8 }_{ -1.4}$| | |$2.7^{+1.3 }_{ -1.0}$| | |$5.5^{+2.6 }_{ -2.0}$| |
| Mw = 2.3 | Mw = 2.2 | Mw = 2.4 | |
|---|---|---|---|
| κ = 0.21 (Brune 1970) | |$21.6^{+10.1 }_{ -7.7}$| | |$15.3^{+7.2 }_{ -5.5}$| | |$30.5^{+14.3 }_{ -10.9}$| |
| κ = 0.26 (Kaneko & Shearer 2014, 2015) | |$11.4^{+5.3 }_{ -4.1}$| | |$8.0^{+3.8 }_{ -2.9}$| | |$16.0^{+7.5 }_{ -5.7}$| |
| κ = 0.28 (Kaneko & Shearer 2014, 2015) | |$9.1^{+4.27 }_{ -3.25}$| | |$6.4^{+3.0 }_{ -2.3}$| | |$12.8^{+6.0 }_{ -4.6}$| |
| κ = 0.32 (Huang et al.2017) | |$6.1^{+2.9 }_{ -2.2}$| | |$4.3^{+2.0 }_{ -1.5}$| | |$8.6^{+4.0 }_{ -3.1}$| |
| κ = 0.372 (Madariaga 1976) | |$3.9^{+1.8 }_{ -1.4}$| | |$2.7^{+1.3 }_{ -1.0}$| | |$5.5^{+2.6 }_{ -2.0}$| |
The stress-drop estimates are derived from |$\Delta \sigma = \frac{7}{16}\left(\frac{\bar{f_c}}{\kappa \beta }\right)^3M_0$| with β as local shearer wave speed 3.26 km s−1 and |$\bar{f_c}$| as 17.1 Hz. κ = 0.26 and κ = 0.28 are suggested for asymmetric elliptical (κ = 0.26) and circular (κ = 0.28) rupture scenarios in Kaneko & Shearer (2014, 2015). The upper and lower stress-drops are estimated with the corner frequency SD (2.3 Hz) in Fig. 11.
Stress-drop estimates (Δσ MPa) from corner frequency.
| Mw = 2.3 | Mw = 2.2 | Mw = 2.4 | |
|---|---|---|---|
| κ = 0.21 (Brune 1970) | |$21.6^{+10.1 }_{ -7.7}$| | |$15.3^{+7.2 }_{ -5.5}$| | |$30.5^{+14.3 }_{ -10.9}$| |
| κ = 0.26 (Kaneko & Shearer 2014, 2015) | |$11.4^{+5.3 }_{ -4.1}$| | |$8.0^{+3.8 }_{ -2.9}$| | |$16.0^{+7.5 }_{ -5.7}$| |
| κ = 0.28 (Kaneko & Shearer 2014, 2015) | |$9.1^{+4.27 }_{ -3.25}$| | |$6.4^{+3.0 }_{ -2.3}$| | |$12.8^{+6.0 }_{ -4.6}$| |
| κ = 0.32 (Huang et al.2017) | |$6.1^{+2.9 }_{ -2.2}$| | |$4.3^{+2.0 }_{ -1.5}$| | |$8.6^{+4.0 }_{ -3.1}$| |
| κ = 0.372 (Madariaga 1976) | |$3.9^{+1.8 }_{ -1.4}$| | |$2.7^{+1.3 }_{ -1.0}$| | |$5.5^{+2.6 }_{ -2.0}$| |
| Mw = 2.3 | Mw = 2.2 | Mw = 2.4 | |
|---|---|---|---|
| κ = 0.21 (Brune 1970) | |$21.6^{+10.1 }_{ -7.7}$| | |$15.3^{+7.2 }_{ -5.5}$| | |$30.5^{+14.3 }_{ -10.9}$| |
| κ = 0.26 (Kaneko & Shearer 2014, 2015) | |$11.4^{+5.3 }_{ -4.1}$| | |$8.0^{+3.8 }_{ -2.9}$| | |$16.0^{+7.5 }_{ -5.7}$| |
| κ = 0.28 (Kaneko & Shearer 2014, 2015) | |$9.1^{+4.27 }_{ -3.25}$| | |$6.4^{+3.0 }_{ -2.3}$| | |$12.8^{+6.0 }_{ -4.6}$| |
| κ = 0.32 (Huang et al.2017) | |$6.1^{+2.9 }_{ -2.2}$| | |$4.3^{+2.0 }_{ -1.5}$| | |$8.6^{+4.0 }_{ -3.1}$| |
| κ = 0.372 (Madariaga 1976) | |$3.9^{+1.8 }_{ -1.4}$| | |$2.7^{+1.3 }_{ -1.0}$| | |$5.5^{+2.6 }_{ -2.0}$| |
The stress-drop estimates are derived from |$\Delta \sigma = \frac{7}{16}\left(\frac{\bar{f_c}}{\kappa \beta }\right)^3M_0$| with β as local shearer wave speed 3.26 km s−1 and |$\bar{f_c}$| as 17.1 Hz. κ = 0.26 and κ = 0.28 are suggested for asymmetric elliptical (κ = 0.26) and circular (κ = 0.28) rupture scenarios in Kaneko & Shearer (2014, 2015). The upper and lower stress-drops are estimated with the corner frequency SD (2.3 Hz) in Fig. 11.
The stress-drop of the M2 earthquake (either from second moments or from corner frequency) is slightly higher but not significantly different from the median of stress-drops of tectonic earthquakes with similar magnitudes at other regions (e.g. Allmann & Shearer 2007; Boyd et al.2017). Similar stress-drop estimates for potentially induced earthquakes have been reported with medians ranging from ∼5 to ∼15 MPa (Clerc et al.2016; Huang et al.2016, 2017; Wang et al.2016; Zhang et al.2016; Boyd et al.2017; Wu et al.2018). The stress-drop estimate in conjunction with the resolved finite source properties do not suggest any abnormal rupture behaviours of this M2 earthquake comparing to tectonic earthquakes. Therefore, this earthquake can serve as a potential example showing induced earthquakes as indistinguishable from natural tectonic earthquakes. On the other hand, the observed stress-drop value here is significantly different than the median values reported in other studies in nearby regions (e.g. Sumy et al.2017; Trugman et al.2017). Sumy et al. (2017) and Trugman et al. (2017) find low stress-drops of small induced earthquakes in Oklahoma and Kansas (∼0.2 to ∼0.4 MPa), which scale with earthquake moments and depths. This difference is unlikely from earthquake location errors, focal mechanism errors, or extreme rupture scenarios (Figs 8 and 9, Table 2). One possible source of the high stress-drop estimate comes from the inaccurate moment magnitude of the earthquake. An Mw 1.3 earthquake would satisfy the resolved optimal source dimension with a 0.23 MPa stress-drop, or an Mw 1.5 earthquake with a 0.46 MPa. However, such a large magnitude deviation for the Mw 2.3 earthquake is unlikely to occur as shown in Sumy et al. (2017). If the local magnitude scales as |${\rm \mathit{ M}_w}=\frac{2}{3}{\rm \mathit{ M}_L}+1$| (Huang et al.2016), the moment magnitude of the M2 earthquake would be Mw 2.53, leading to a stress-drop of 16.2 MPa correspondingly. Possible physical mechanisms explain such variability often involves fluids, for example, fluid pressure at the interface reducing the normal stress, and stable to unstable slip transition on given asperities due to presents of fluids (e.g. Lengliné et al.2014). However, as shown by reservoir induced earthquakes, the presence of water do not affect stress-drop (e.g. Tomic et al.2009). With only one event in this study, it is difficult to conclude whether the observed stress-drop is representative for earthquakes in the region. The unanswered question warrants continued monitoring and further investigation of earthquakes in the region.
To further explore the observed spatial patterns of apparent duration, S-wave corner frequency and fall-off exponent, we compare the observations with a synthetic asymmetric elliptical source scenario evaluated in McGuire & Kaneko (2018; Fig. 12). The simulated earthquake mostly unilaterally propagated at 90 per cent of local shear wave speed with Mw 3.9 (Fig. 12). The apparent durations of the simulated earthquake are smooth and vary gradually with incident angles, while the corner frequencies and fall-off exponents show great variabilities over the focal sphere. The corner frequencies are larger for receivers in the rupture direction than the region away from rupture propagation, while fall-off exponents show complex spatial pattern (Figs 12D and F). Although the simulated earthquake magnitude is different than that of the M2 earthquake, apparent duration, S-wave corner frequency and fall-off exponent share similar spatial patterns between the model and the real earthquake. These general agreements validate the dynamic rupture simulation, and highlight that rupture directivity can be well resolved with dense nodal array even for subshear small earthquakes.
Comparison between the observed (left column) and synthetic (right column) apparent duration, S-wave corner frequency and fall-off exponents. The synthetic simulation is from Kaneko & Shearer (2014, 2015). The simulated earthquake is an asymmetric elliptical source that unilaterally propagated at 90 per cent of local shear wave speed with Mw 3.9. Note that the observed and simulated earthquake magnitudes are different.
Seismic instrumentation has been evolving along with the forefront of seismology research. Large-N nodal arrays suggest one future direction for data acquisition, which are composed of hundreds to thousands short period seismic sensors deployed with tens of kilometres (e.g. Schmandt & Clayton 2013; Sweet et al.2018). These large-N nodal arrays have been implemented to study local seismic structures (e.g. Schmandt & Clayton 2013; Nakata & Beroza 2015; Ward & Lin 2017), infer fault zone geology (e.g. Qin et al.2018), detect seismic events (e.g. Inbal et al.2016; Li et al.2018; Meng & Ben-Zion 2018) and track anthropogenic footprints (e.g. Riahi & Gerstoft 2015). Here, we show that second moment method can be implemented with the large-N nodal arrays to resolve small earthquake rupture dynamics at great precision. The rupture area derived from the second moment method can be used to estimate small earthquake stress-drop directly, which is independent of any a priori model. The nodal systems of the IRIS wavefield experiment have analogue-to-digital converter resolution as 24 bits, preamplifier gain adjustable from 0 to 36 dB, and sampling rate adjustable from 250 to 2000 Hz. These flexibilities of the nodal instrument suggest the system can be deployed to observe earthquakes with a wide range of magnitude up to M3 earthquakes. Deployment of such dense nodal arrays on active faults can potentially aid in understanding earthquake scaling and help in understanding the large scatter of earthquake stress-drop estimates.
6 CONCLUSIONS
The IRIS Community Wavefield Experiment in Oklahoma provided a great opportunity to examine local small earthquake finite source properties in high resolution. In this study, the nodal array data are implemented to resolve the finite source attributes of one local M2 earthquake by the second moment method. In total, seven EGFs and 658 SH apparent duration measurements were used to invert for the second moments. The M2 earthquake unilaterally ruptured with Lc ∼70 m and Wc ∼45 m at 90 per cent local S-wave speed (2.93 km s−1), lasting for about 0.019 s. With the derived rupture area, we obtain the stress-drop of the earthquake as 7.3 MPa, which agrees with stress-drop estimates of other similar magnitude earthquakes hosted by natural faults. The kinematic and dynamic properties of this M2 earthquake suggest little difference in its rupture process from typical large earthquakes. In addition, we fit the Boatwright spectral model to investigate the spatial variation of the corner frequency and fall-off exponent of the earthquake. The observed spectral parameters show clear directivity effects of the rupture propagation. By applying the second moment method to nodal array data, we can obtain robust estimates of key finite fault parameters, which can help in understanding small earthquake dynamics and may have potential to address the earthquake scaling problem. The tight constraints on rupture area from this array are noteworthy in that the array geometry (three short lines) was not ideal for covering the focal sphere even for events located in its centre. Future deployments of three-component nodal arrays, perhaps with an even larger number of instruments and with a more appropriate geometry, could improve on constraining rupture area for small to moderate earthquakes.
SUPPORTING INFORMATION
Supplementary data are available at GJI online.
Figure S1. Vertical component of the three seismic lines (P wave). The records are filtered at 10–100 Hz. The red crosses indicate the starting point of a time window derived from cross-correlation.
Figure S2. Transverse component of the three seismic lines (SH wave). The legends are similar to Fig. S1.
Figure S3. Radial component of the three seismic lines (Sv wave). The legends are similar to Fig. S1.
Figure S4. Seven EGFs of Line 1 (SH wave). The redlines indicate the time window used for cross-correlation with the main-shock records. The EGF colours correspond to the measurement colour scheme in Fig. 4.
Figure S5. Seven EGFs of Lines 2 and 3 (SH wave). The legends are similar to Fig. S4.
Figure S6. Distributions of characteristic length (A), width (B), duration (C) and rupture velocity (D). The distributions are derived from 5000 times bootstrap resampling including the data outliers shown in Fig. 4(B).
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.
ACKNOWLEDGEMENTS
We thank the editor Martin Mai and the two reviewers for their constructive suggestions, which have led to improvements in the paper. We thank Yoshihiro Kaneko for sharing his dynamic rupture models, Xiaowei Chen, Heather DeShon and Louis Quinones for sharing their earthquake catalogues, Maria Beatrice Magnani for sharing local velocity model, and Yihe Huang and Jianhua Gong for helpful discussions. The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata and/or derived products used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681. WF is currently supported by the Postdoctoral Scholar Program at the Woods Hole Oceanographic Institution, with funding provided by the Weston Howland Jr. Postdoctoral Scholarship. JM was partially supported by SCEC grant #17177 at Woods Hole Oceanographic Institution. This research was supported by the Southern California Earthquake Center (Contribution No. 8014). SCEC is funded by NSF Cooperative Agreement EAR-1033462 and USGS Cooperative Agreement G12AC20038. All the processed data are available upon request to the authors.
