Automatic determination of focal depth with the optimal period of Rayleigh wave amplitude spectra at local distances

SUMMARY

Focal depth of earthquakes is essential for studies of seismogenic processes and seismic hazards. In regions with dense seismic networks, focal depth can be resolved precisely based on the traveltime of P and S, which is less feasible in case of sparse networks. Instead, surface waves are usually the strongest seismic phases at local and regional distances, and its excitation is sensitive to source depth, thus theoretically important for estimating focal depth even with a limited number of seismic stations. In this study, short-period (0.5–20 s) Rayleigh waves are explored to constrain focal depths. We observe that the optimal period (the period corresponding to the maximum amplitude) of Rayleigh waves at local distances (≤200 km) shows an almost linear correlation with focal depth. Based on this finding, we propose an automated method for resolving the focal depth of local earthquakes using the linear regression relation between the optimal period of Rayleigh wave amplitude spectra and focal depth. Synthetic tests indicate the robustness of this method against source parameters (focal mechanism, source duration and non-double-couple component) and crustal velocity structure. Although the attenuation (Q factor) of shallow crust can introduce complexities in determining focal depth, it can be simultaneously estimated if a sufficient number of stations are available. The proposed method is applied to tens of small-to-moderate earthquakes (M_w 3.5–5.0) in diverse tectonic settings, including locations in the United States (Oklahoma, South Carolina, California, Utah, etc.) and China (Sichuan, Shandong). Results demonstrate that reliable focal depth, with uncertainty of 1–2 km, can be determined even with one or a few seismic stations. This highlights the applicability of the method in scenarios characterized by sparse network coverage or historical events.

1 INTRODUCTION

Focal depth of earthquakes provides key information to reveal the seismogenic structure (i.e. the depth extent of active faults, e.g. Shearer 1997; Ross et al. 2017) and crustal rheology structure (i.e. brittle-ductile transition, e.g. Scholz 2019), infer the rupture nucleation or co-seismic rupture propagation process (e.g. Dodge et al. 1996; Davenport et al. 2015) and evaluate the seismic hazard (e.g. Papazachos et al. 1993). Besides, focal depth is a crucial criterion in identifying human-induced earthquakes, especially for regions where historical seismicity is low or tectonic earthquakes occur at depth ∼10 km (e.g. Davis & Frohlich 1993; Keranen et al. 2013; Chen et al. 2017).

Available methods to resolve the focal depths can be generally classified into two groups, which are based on the traveltime and seismic waveforms, respectively. The conventional earthquake location methods based on the absolute or differential traveltime of Pg and Sg can well constrain the focal depth in regions where there is at least one station with the epicentral distance of less than twice the focal depth and the reliable 3-D velocity structure is available (Thurber et al. 2006). However, in the case of sparse seismic networks, the travel time or the propagation distance of direct waves is less sensitive to focal depth and these methods might not work well. Alternative approaches, such as the regional seismic depth phases, regional waveform inversion and surface wave amplitude features, make full use of the abundant information of seismic waveforms and help to better constrain the focal depth.

The seismic depth phases reflected from the free surface provide effective constraints on focal depth. For an epicentral distance of 30–300 km, focal depth can be estimated by modelling depth phases sPL, sPg, sPmP and sPn, etc. (Ma & Atkinson 2006; Ma 2010; He et al. 2019). Yuan et al. (2020) combined possible depth phases and scan waveforms that match depth-phase templates, which improves the efficiency and accuracy of the depth‐phase methods. For moderate earthquakes, the focal depth can be estimated together with other source parameters (e.g. focal mechanism) via regional waveform inversion (Cesca et al. 2010), while for small earthquakes, amplitude calibration or a 3-D velocity model is required to model short-period waveforms (Saikia et al. 2001; Wang & Zhan 2020).

Usually, surface waves are the strongest seismic phases at local and teleseismic distances, and the excitation is sensitive to source depth due to the approximately exponential decay of surface wave displacements versus depth. Thus, surface wave amplitude features (e.g. spectral nulls, the amplitude ratio between surface and body waves, as well as the whole amplitude spectra) can be used to estimate the focal depth (Tsai & Aki 1971; Nguyen & Herrmann 1992). As some of these features might also be affected by the focal mechanism, the focal depth is usually resolved together with the focal mechanism. Fox et al. (2012) inverted the earthquake source parameters by fitting the intermediate-period (15–100 s) surface wave amplitude spectra at regional and teleseismic stations, and determined the focal mechanism, depth and scalar seismic moment of 20 M_w 4.3–6.4 earthquakes. Heyburn et al. (2013) combined the surface wave amplitude spectra and teleseismic depth phase (pP) to estimate the source depth. Jia et al. (2017) jointly inverted the source depth of M_w 5.6–6.0 earthquakes with local seismic waveforms and regional Rayleigh wave amplitude spectra (10–100 s). In these studies, amplitude spectra of intermediate-period (10–100 s) surface waves at regional distances are utilized, while surface waves at shorter periods have not been well explored yet, which can provide valuable information on determining focal depths. Indeed, studies of short-period surface waves at regional distances (∼1000 km) are challenging due to either strong scattering or substantial attenuation, which lead to quite weak signals. But at local distances (≤200 km), short-period surface waves are expected to be well observable for crustal earthquakes. For example, in recent years, multiple large-scale broad-band seismic arrays with station spacing of 50–100 km have been deployed world-widely, such as PASSCAL, USArray and ChinArray, providing abundant local waveform data and making it possible to explore the potential of short-period surface waves in constraining focal depth. For instance, short-period Rayleigh waves are observed for a shallow event (depth=3 km) in Sichuan at an epicentral distance of 15 km (Luo et al. 2010), but the application of short-period Rayleigh waves to quantitatively constrain focal depth requires more investigations and validations.

In this paper, we propose a novel automated method for resolving the focal depth of earthquakes using the optimal period, which corresponds to the maximum of the short-period (0.5–20 s) Rayleigh wave amplitude spectra (referred to as optimal period). To assess the effectiveness of the proposed method, extensive synthetic tests are conducted, evaluating its robustness against various source parameters such as focal mechanism, source duration and non-double-couple components, as well as crustal structures including velocity and attenuation. Furthermore, the focal depths of tens of earthquakes (M_w 3.5–5.0) in different tectonic settings across the United States (Oklahoma, South Carolina, California, Utah, etc.) and China (Sichuan, Shandong) are determined, and the applicability of this method is thoroughly investigated within these diverse tectonic environments.

2 METHODOLOGY

2.1. Optimal period of Rayleigh wave amplitude spectra

For a double-couple seismic source at a depth of h, the excited vertical displacement of Rayleigh waves propagating in a 1-D Earth or ‘layered model’ can be written as below (Aki & Richards 2002).

where

where M_ij (i,j = x, y and z) are moment tensor elements describing the fault geometry and earthquake magnitude, r represents the epicentral distance between the source and a station, and |$\phi $| is the azimuth. The phase velocity c and group velocity U, which are functions of wavenumber k or the wave period, together characterize how fast surface waves propagate along the surface. Note that they depend on only the Earth model (i.e. density ρ and seismic velocity) and has nothing to do with source property. Similarly, the eigenfunctions r₁ and r₂ are also functions of wave period and their shapes are independent on the moment tensor. They describe the amplitude and phase information of surface wave across the depth and play a key role in our proposed method. The focal depth appears only in the eigenfunctions or their depth derivatives in eq. 1, so the eigenfunctions act as ‘kernels’ to link the focal depth h and displacement u_z. With eq. (1), one can infer the focal depth together with other source parameters (i.e. moment tensor and horizontal locations) by fitting the observed displacement waveforms at multiple periods. However, the waveforms are not only dependent on the source but also affected by the accumulated effects of wave propagation. These effects cannot be accurately predicted without a well-constrained shallow Earth structures, which is the common case for many regions, and thus usually become the largest uncertainty source of the depth estimation.

The amplitude spectra pattern in the short periods (e.g. 0–20 s) are most sensitive to the depth of crustal earthquakes. This high sensitivity is primarily due to the large gradient of eigenfunctions with the depth (Fig. S1, Supporting Information). In this study, we propose to use the amplitude spectra across different periods to estimate the focal depth. The idea of this method is to extract the depth-sensitive amplitude spectra from observed waveforms for a depth estimation while mitigate other ‘contaminating’ effects in the waveforms due to other source parameters or Earth model uncertainty. For instance, the phase information of waveforms is prone to be affected by the Earth model uncertainty, the absolute amplitude is predominantly influenced by the earthquake magnitude, and the spectral null is sensitive to the focal mechanism. In contrast with these depth-insensitive or error-prone variables, we find that the period corresponding to the maximum amplitude, referred to as ‘optimal period’, serves as a robust indicator of focal depth.

To demonstrate this relationship between the optimal period and the focal depth, we make a four-layer model (sedimentary layer, upper crust, lower crust and the mantle half-space, Table S1, Supporting Information) and use the CPS (Computer Programs in Seismology) package (Herrmann 2013) to compute the fundamental mode spectra for a depth range of 0–20 km. Here we show the amplitude spectra that are each independently normalized by their maxima (Fig. 1). It is clear that the amplitude spectra pattern democratically changes with the depth for all the typical focal mechanisms, which can be explained by the large gradient of eigenfunctions with depth (Fig. S1, Supporting Information). However, the Rayleigh wave amplitude spectra pattern also strongly depend on the source mechanism. For example, the spectral null for strike-slip faulting and thrust faulting can be observed and the period of the null increases as the earthquake becomes deeper. However, for vertical dip-slip faulting, the spectral null cannot be observed at this azimuth. Compared with the spectral null, a more robust feature is that the period with the largest amplitude or energy (referred to as the optimal period) increases stably as the source becomes deeper. We measure the optimal period from the spectra via fitting the spectra near the maxima with a parabola function, which provides a more precise estimation of the location of the maximum point and is more robust to outliers (e.g. Chen et al. 2015). We find that the optimal period is almost linearly correlated to the focal depth for all the cases (Figs 1d–f), thus offering a helpful tool for determining focal depth. For a given reference velocity model and focal mechanism, the nearly linear relationship between the optimal period and focal depth can be constructed, then the focal depth can be resolved once the optimal period is measured (further details regarding this process can be found in Section 3). It should be noted that in certain scenarios, such as deeper events or the azimuths near the nodal, the measurement of the optimal period might experience more uncertainty. This can be attributed to the small absolute amplitude (Fig. S2, Supporting Information), which makes it susceptible to background noise. In such cases, it would be better to exclude waveforms with a low signal-to-noise ratio (SNR) when applying the method to real data.

2.2 Synthetic tests based on eigenfunctions

To investigate the robustness of this method, we conduct a series of synthetic tests for different source parameters (focal mechanism, source duration and non-double-couple component) and crustal structures (velocity and attenuation).

First, we investigate the influence of focal mechanism and azimuth on the optimal period. To simulate more realistic cases, we select the focal mechanisms from four events (M_w 4.0–4.5) in the United States with different tectonic settings (Oklahoma, South Carolina, California and Utah) (Fig. 2 and Table 1), including strike-slip faulting, thrust faulting and normal faulting. The CUS velocity model (Table S2, Supporting Information) is used as the 1-D reference velocity model, and the quality factors of P wave (Qp) and S wave (Qs) are 1000 and 500, respectively.

The optimal periods for different focal mechanisms and azimuths are calculated for depths of 5, 10 and 15 km. For the case of strike-slip faulting, the azimuthal variation of the optimal period is minor and negligible (Fig. 3a). For the case of dip-slip faulting, both the mean value and the azimuthal variation are larger than that of strike-slip faulting (Figs 3b–d). The variations of optimal periods are within ± 10 per cent of the mean value, implying the maximum uncertainty of 10 per cent in the resolved depth, that is, 1 km for source at depth of 10 km. For real cases, this uncertainty can be suppressed by averaging the optimal methods at stations of different azimuths.

For real cases, it is possible that the focal mechanisms are not perfectly known, so we conduct forward tests to quantify the errors that may arise due to inaccurate focal mechanisms. The uncertainty in focal mechanism solutions derived with different methods or from different catalogues is usually smaller than 15° (Helffrich 1997; Frohlich & Davis 1999), but occasionally, the difference can be up to 30° (Duputel et al. 2012). Thus, for each of the four events, we generate 12 groups of different focal mechanisms, which are respectively 15° and 30° variation in strike, dip and rake for NP1, to calculate the theoretical optimal period and investigate the sensitivities on the uncertainties in focal mechanism solutions.

The results (Fig. S3, Supporting Information) show that (1) the distribution of the bias is not symmetrical with respect to that of focal mechanism, and the percentage of the bias does not change significantly with the focal depth; (2) for a variation of 15°, the bias of 83 per cent cases (397 in 480 cases) are within 5 and 98 per cent cases (472 in 480 cases) are within 10 per cent, and all the cases are within 15 per cent; (3) for a variation of 30°, the bias of the 73 per cent cases (351 in 480 cases) are within 10 and 89 per cent cases (428 in 480 cases) are within 15 per cent, and all the cases are within 25 per cent. Overall, for most cases, a variation of 15° in focal mechanism solution would result in less than 10 per cent bias in the optimal period, which propagate into 10 per cent bias of the resolved focal depth.

Next, we investigate the sensitivity of this method to the crustal velocity model. Taking the strike-slip faulting mechanism as an example, we adopt the CUS model as the reference velocity model and construct 100 models via varying the CUS model by adding random perturbations (≤5 per cent) (Fig. 4a). The results show that for input depth range of 1–17 km, the bias in estimated depth is smaller than 2 km for 95 per cent of cases, and for a deeper source (18–20 km), the bias is smaller than 2 km for 70 per cent of cases (Fig. 4b), indicating the method is relatively robustness when the velocity model is relatively well constrained. Then, we carry out tests with two models that are distinct enough, including a basin-range crustal model (WUS) and a stable shield model (CUS). The difference between the two velocity models is largest near the ground surface, which is 1.6 km s⁻¹ (47 per cent) for P wave and 0.9 km s⁻¹ (45 per cent) for S wave (Fig. S4a, Supporting Information). We calculate the eigenfunctions and find the difference between the calculated optimal period for two velocity models ranges from −0.8 to 0.9 s (Fig. S4b, Supporting Information). The largest relative difference occurs at 1 km (−54 per cent), and the relative difference for depth of 5–20 km is within ± 10 per cent. For all the cases, the difference in the estimated focal depth is smaller than 2 km and the average difference is 0.9 km, suggesting that a large uncertainty in crustal velocity model may result in 1–2 km bias of the resolved focal depth.

Short-period Rayleigh waves usually attenuate for long propagation distance, and the optimal period may be shifted towards longer period for regions with high attenuation. Thus, we carry out sensitive tests on the Q factor and evaluate the influence of attenuation of shallow crusts. Three tests are performed, in which the range of Qp is 10–500 (Fig. 5). In Case 1, the Q factor for depth shallower than 5 km varies from 20 to 500. Here, the focal depth is set to be 5 km. The results show that if Qp is larger than 140, the optimal period variation for distances ranging from 20 to 100 km is smaller than 0.5 s. For high attenuation (Qp < 40), the variation is larger than 1.0 s. Even for the nearest station, the attenuation leads to an increase of 0.5–1.0 s in the optimal period. For Case 2, the source depth is 10 km. For Case 3, the source depth is 10 km and the Q factor for depth shallower than 10 km is tested. For these three cases, we calculate the slope of the optimal period variation (change of optimal period due to longer distances, with dimension in s/km) and find that the slopes of Cases 1 and 3 are similar, while the slope of Case 2 is smaller. This sensitivity test suggests that the averaged Q factor of shallow crust can be estimated roughly with the variation of the optimal period, if sufficient number of seismic stations are available.

2.3 Synthetic tests based on waveforms

Previous synthetic tests are based on eigenfunctions, which are calculated directly for the given source parameters and crustal velocity models. To verify the above results, we also calculate synthetic seismograms with the F–K method (Zhu & Rivera 2002) and measure the optimal period from the waveforms. As an example, we calculate the displacement synthetics with the focal mechanism of the 2014 Oklahoma event (10°/86°/−180°) for depth of 1–20 km and distance of 30–100 km. For depth of 5 km, Rayleigh waves are clear and observable for distance ranging from 30 to 100 km. From visual inspection, the optimal periods for different distances are similar (Fig. 6). We apply the do_mft algorithm in the CPS package (Herrmann 2013) to the waveforms and get the energy or amplitude for different periods. Then, we fit the spectra near maxima with a parabola function and determine the optimal period. The measured optimal periods for different distances are consistent with each other and are also consistent with the calculated optimal period based on eigenfunctions (Figs 6c and d).

We measure the optimal periods for different depths and distances, and calculate the difference between the values calculated with eigenfunctions. The results show that the measured optimal periods for different depths are consistent with that calculated with eigenfunctions (Fig. 7), suggesting that the relationship between the optimal period and source depth can be efficiently calculated based on eigenfunctions. Furthermore, for shallow events, the optimal period can be measured for a large range of epicentral distances, while for deeper events, the optimal period can be measured only at stations tens of kilometers away. This is reasonable as Rayleigh waves develop at distances of about three wavelengths and 3–5 times of focal depth (Aki & Richards 2002). Taking the distance of 75 km as an example, the optimal period does not change for events deeper than 15 km events, indicating the failure for resolving the depth of deeper events at this distance. We find that 1/5 distance is a proper criterion for successful application to observed seismic data. That is, for depth shallower than 1/5 of the distance, the optimal period can be stably measured with a bias smaller than 1 s. This seems an ad-hoc criterion as focal depth is not known yet as soon as an earthquake just occur. This is indeed a problem for earthquakes near subduction zones, where earthquakes can be as deep as a few hundred kilometres. But for crustal earthquakes, seismic events are usually shallower than 30 km, thus epicentral distance around 150 km would be able to resolve the focal depths. As an approximation, for an earthquake with the source depth unresolved, whether the short-period Rayleigh can be observed at close stations provides a preliminary estimation of the focal depth.

Next, we test the validity of this method against more source parameters, including the source duration and non-double-couple components. In this method, reliable measurement of the optimal period of Rayleigh waves is essential, but it could be affected by source duration. For example, the optimal period is ∼3 s for depth of 5 km for an event with very short duration. For an M7 event, the source duration is ∼10 s, which of course would lead to the failure of this method. Here we only test for small-to-moderate earthquakes, including the source duration of 0.3 s (∼M_w 4), 1.0 s (∼M_w 5.0) and 1.7 s (∼M_w 5.5). Synthetics show that the high-frequency phases become weaker when the source duration increases (Fig. 8a). Following similar procedures, we measure the optimal period for different cases. The results show that the difference of most cases is within 0.5 s, but the method would fail for depth of 1–2 km when the duration is 1.7 s (Fig. 8b). Thus, this method can be applied to most small-to-moderate (M_w 3.0–5.5) events, but may fail for extremely shallow moderate earthquakes. This test also implies that this method is relatively robust against the rupture directivity, which would cause the azimuthal variation of spectra. For an M_w 5 event, we assume that the rupture length is 3.0 km and the Rayleigh wave speed is about 3.0 km s⁻¹, then the azimuthal variation of apparent source duration is about 1 s, and would not substantially affect the application of this method.

Volumetric change or fault complexity in the source excitation may lead to non-double-couple components in the source solution, thus leading to some effects on estimating focal depth with the optimal period. We take the Oklahoma event as an example to show the effect of the non-double-couple component. The USGS regional moment tensor solution contains 87 per cent of double couple components (Fig. S5a, Supporting Information), and we calculate the synthetics with the moment tensor for depths of 1–20 km. We find that the bias of optimal period in most cases is within 0.2 s (Fig. S5b, Supporting Information). This test suggests that this method is relatively robust against weak non-double-couple components, which is the case for most small-to-moderate earthquakes.

3 PROCEDURES FOR DETERMINING THE FOCAL DEPTH

Fig. 9 illustrates the flowchart used to estimate focal depth using the optimal period. This process involves two main steps. First, the optimal period is determined by analysing recorded waveforms. Secondly, a theoretical regression relation between the optimal period and focal depth is calculated, and the focal depth is determined using the corresponding regression relations.

In the first step, vertical components of broad-band seismograms are collected at local distances ranging from 20 to 150 km. Next, the instrument responses are deconvolved to obtain the displacement seismogram, followed by the removal of the mean value and linear trend. Manual measurements of the amplitude spectra can be performed using the do_mft algorithm (Herrmann 2013), while automated measurements can be achieved using either the sacmft96 algorithm (Herrmann 2013) or the frequency–time analysis (FTAN) algorithm (Bensen et al. 2007). Then, the optimal period is determined by fitting the spectra near maxima with a parabola function.

In regions with dense networks, the optimal periods at different stations are linearly fitted against the epicentral distance to mitigate the influence of attenuation (eq. 2). The predicted optimal period from the nearest station is utilized as the final optimal period for the earthquake.

where |$dis{t}_i$| and |$o{p}_i$| are the epicentral distance and the measured optimal period at stationi, slope and b can be obtained with the least square fit. |$op\_fit$|is the predicted optimal period at the closest station with distance of |$dis{t}_{\min }$|⁠.

In regions where network coverage is limited, the final optimal period for the earthquake is determined by taking the average of the obtained optimal periods (eq. 3), which ensures a more reliable estimation.

where|$o{p}_i$|is the measured optimal period at station i and|$op\_ave$|is the average optimal period of the measurement at N stations.

In the second step, the process begins by calculating the theoretical regression relation between the optimal period and focal depth. This is achieved by utilizing the focal mechanism and a local crustal velocity model, which can be obtained from global models such as Crust1.0 or Crust2.0 (Laske et al. 2004; Laske et al. 2013), or from previous studies involving receiver functions and other geophysical methods. By combining the crustal velocity model with the focal mechanism of the target event, theoretical amplitude spectra are computed based on eigenfunctions for various azimuths (ranging from 0° to 360° with an interval of 30°) and depths (ranging from 1 to 20 km with an interval of 1 km). From these spectra, the optimal periods are extracted and averaged across different azimuths to obtain a mean value for each depth. Subsequently, a linear least-square fitting is performed to establish the theoretical regression relation between the optimal period and focal depth. This regression relation, along with the optimal period determined in the first step, is used to determine the focal depth of the target earthquake.

In cases where the focal mechanism of the target earthquake is unavailable, historical events or standard faulting mechanisms such as strike-slip faulting (0°/90°/90°) and thrust faulting (0°/45°/90°) can be used as substitutes. To assess the applicability of the regression relation to events with different faulting mechanisms, forward tests are conducted for the most extreme scenario. The theoretical optimal periods are calculated for depths ranging from 0.5 to 20 km using the vertical strike-slip mechanism, then the corresponding depths are derived using the regression relation that is calculated for thrust faulting. The results, as shown in Fig. S6 (Supporting Information), demonstrate that the absolute bias in depth is less than 2.5 km for all cases, and the relative bias is below 20 per cent for depths exceeding 4.5 km.

4 APPLICATION TO OBSERVED EARTHQUAKES

To investigate the applicability of the new method in different tectonic settings, we choose tens of earthquakes (M_w 3.5–5.0) in the United States (Oklahoma, South Carolina, California, Utah, etc.) and China (Sichuan, Shandong) as examples and determine the focal depths of these events. Most of these earthquakes have been well studied with independent methods, providing nice benchmarks for our method.

4.1 The 2014 Oklahoma M_w 4.2 event

Oklahoma has experienced a sharp increase in seismicity since 2009 and induced earthquakes have become an important topic (e.g. Ellsworth 2013). The 2014 October 10 M_w 4.2 earthquake is a strike-slip event in the shallow crust (Table 1), with source parameters well studied by McNamara et al. (2015) and Yuan et al. (2020) using data sets from local portable stations and regional seismic networks.

We collect the broad-band waveforms at 32 local stations with epicentral distances of 20–100 km, and retain waveform data from 31 stations with high SNR. We check the waveforms visually and find the group velocity of Rayleigh waves is smaller than 2.5 km s⁻¹ (Fig. 10a). Polarization analysis with Hilbert transform on the radial component confirms that Rayleigh waves are correctly identified. The amplitude spectra of all the stations show high consistency (Fig. 10c). Considering the realistic error in amplitude measurements, we fit a parabola function to five points near the maxima, instead of three, to measure the optimal period. The five points (instead of three points) near the maxima with a parabola function to measure the optimal period. Furthermore, we estimate the error associated with this measurement by providing the standard deviation (STD) obtained from the fitting process. The optimal periods for most stations are within 2.0–3.0 s with the mean value of 2.53 s, and the averaged STD is 0.022 s (Fig. S7, Supporting Information). The optimal periods increase slightly for larger distances, thus we fit the data linearly and the predicted optimal period at the closest station is 2.32 s. The fitting slope is 0.0051 s km⁻¹, corresponding to Qp of ∼200 (Fig. 5), which is similar to the results from the crustal seismic (Lg) attenuation tomography (Levandowski et al. 2021).

Here, we perform forward tests on how the averaging scheme works. We randomly select a subset of stations from available stations and calculate the average optimal period. Each test on the number of the selected stations is repeated for 100 times (Fig. S8a, Supporting Information). We calculate the STD of each case, and adopt the difference between the average optimal period for all the stations and the optimal period predicted at the closest station as a criterion or reference to assess how the averaging scheme works (Fig. S8b, Supporting Information). The result shows that the STD is smaller than the reference if there are at least three stations, indicating the applicability in sparse networks.

Following the procedures in Fig. 9(b), we calculate the theoretical optimal periods for depth 1–20 km and azimuth 0°–360° with the CUS model. Then, we average the optimal periods for different azimuths and fit the data linearly against depth to get the coefficients. The slope is 1.759 km s⁻¹ and the intercept is −0.065 km. We combine the fitting coefficients and the observed optimal period (2.32 s for the best-fitting optimal period and 2.53 s for the averaged optimal period) to determine the depth. The resolved depth is 4.0 km for the best-fitting optimal period and 4.4 km for the averaged optimal period (Table 2). The depth is consistent with the values obtained with the depth-phases methods (5.1 ± 0.8 km, Yuan et al. 2020) and waveform inversion (4 km, t. Louis University Earthquake Center solution, referred to as SLU solution).

4.2 The 2014 M_w 4.1 South Carolina event

The 2014 February 15 M_w 4.1 South Carolina earthquake is a shallow thrust event near the border between Georgia and South Carolina (Fig. 2 and Table 1), and its focal depth has been studied by modelling the local depth phases (Daniels et al. 2020; Yuan et al. 2020).

We collect the broad-band waveforms at 23 local stations and retain waveform data from 22 stations with high SNR. The optimal periods for most stations show a significant linear increase against epicentral distance (Fig. 11). The predicted optimal period at the closest station is 2.64 s, while the mean optimal period is 3.30 s and the averaged STD is 0.003 s. In the linear fitting, the slope is 0.0213 s km⁻¹, implying high attenuation. Although this feature has not been found in Lg attenuation tomography (Levandowski et al. 2021), the high attenuation in Atlantic Coastal Plain sediments has been suggested from reflections and sediment-guided P waves (Chapman et al. 2008; Guo & Chapman 2019). Then, we calculate the theoretical regression relation with the CUS model. The slope is 1.458 km s⁻¹ and the intercept is 0.161 km. The resolved depth is 4.0 km for the best-fitting optimal period and 5.0 km for the averaged optimal period (Table 2). The shallower depth (4.0 km) would be more reliable due to the high attenuation. The depth is consistent with depth-phases methods (3.7 ± 0.7 km, Yuan et al. 2020; 4 km, Daniels et al. 2020) and waveform inversion (5 km, SLU catalogue).

4.3 The 2020 M_w 4.5 California event

We also apply this method to moderate earthquakes in the middle crust. The 2020 September 19 M_w 4.5 California earthquake is a thrust event in the Los Angeles Basin where the seismic network is dense (Fig. 2 and Table 1), and the focal depth has been well constrained with waveform inversion demonstrating that the earthquake occurred at depth about 16 km (USGS regional moment tensor solution).

For this event, the minimum epicentral distance that Rayleigh waves can be observed is larger than shallower events, thus we collect the broad-band waveforms at stations with epicentral distances of 60–110 km. The optimal periods for most stations show a significant linear increase against epicentral distance, and are more scattered around distance of ∼60 km (Fig. 12), which may be caused by the less developed surface waves at short distances. The predicted optimal period at the closest station is 10.20 s. The mean optimal period is 11.38 s and the averaged STD is 0.001 s. The slope is 0.042 s km⁻¹, corresponding to strong anelastic attenuation, which is also noted in previous studies (e.g. Olsen et al. 2003; Liu et al. 2021). Then, we calculate the theoretical regression relation with the standard 1-D Southern California velocity model (SoCal, Hadley & Kanamori 1977). The slope is 1.567 km s⁻¹ and the intercept is −0.019 km. Then, the resolved depth is 15.9 km for the best-fitting optimal period and 17.8 km for the averaged optimal period (Table 2). Similar to the South Carolina earthquake, the shallower depth (15.9 km) is more reliable, which is consistent with the results from waveform inversion (16 km, USGS regional moment tensor solution).

4.4 The 2020 M_w 4.0 Utah event

The 2020 April 15 M_w 4.0 Utah earthquake is a normal event in the upper crust and is located where the seismic network is relatively sparse (Fig. 2 and Table 1). We collect the broad-band waveforms at 14 local stations. The linear fitting, the predicted optimal period at the closest station is 5.93 s, while the mean optimal period is 6.48 s and the averaged STD is 0.002 s (Fig. 13). The slope (0.018 s km⁻¹) is less constrained due to the sparse sampling data. Following similar procedures, we calculate the theoretical regression relation with the western United States (WUS) model (Herrmann et al. 2011). The slope is 1.637 km s⁻¹ and the intercept is −0.395 km. The resolved depth is 9.3 km for the best-fitting optimal period and 10.2 km for the averaged optimal period (Table 2). The depth is consistent with USGS regional moment tensor solution (9 km) and SLU solution (11 km).

4.5 The 2019 M_w 4.0–4.3 Weiyuan sequence and the 2022 M_w 3.7 Weifang earthquake

We also apply this method to several small to moderate earthquakes in China. The 2019 Weiyuan earthquake sequence in Sichuan Basin occurs at extremely shallow depths and has been inferred to be related to injection activities (Yang et al. 2020). This sequence includes three M_w 4.0–4.2 events (Table S3, Supporting Information). There are three broad-band stations within 100 km, and the amplitude spectra of the stations show high consistency (Fig. S9, Supporting Information). The 1-D velocity model from Zhao et al. (1997) is used as the reference model, which works well in source inversion (Yang et al. 2020; Yi et al. 2020). Focal mechanisms from waveform inversion in Yang et al. (2020) are used for this study (Table S3, Supporting Information). Since there are only three stations, the average optimal periods for three events (1.99, 1.97 and 1.74 s) are used to determine the depth. The resolved depths for three events are 3.0, 3.0 and 2.7 km, respectively, which is consistent with the depth determined with waveform inversion (2.5–2.6, 2.5–2.9 and 1.0–2.0 km, Yang et al. 2020; Yi et al. 2020).

The 2022 May 1 M_w 3.7 Weifang earthquake is a normal event in Shandong province. We collect the broad-band waveforms at 17 local stations (Fig. S10a, Supporting Information). The predicted optimal period at the closest station is 1.99 s and the mean optimal period is 2.35 s. Adopting the Crust2.0 model (Laske et al. 2004) as the reference velocity model, the resolved focal depth is 2.6 km for the best-fitting optimal period and 3.1 km for the averaged optimal period. We also invert the source parameters with the Cut-And-Paste method (Zhu & Helmberger 1996) and broad-band waveforms at 23 local stations. The focal depth obtained in the waveform inversion is about 3 km (Fig. S10d, Supporting Information), consistent with the optimal period method.

5 DISCUSSION

To further investigate the stability and reliability of the method, here we apply this method to M3.5 + earthquakes in the CEUS. There are 21 M3.5 + earthquakes during 2014–2021 with source parameters resolved in the SLU catalogue (Fig. 14a and Table S5, Supporting Information). Similar procedures (Fig. 9a) are adopted to get the observed optimal period and focal depth. In the forward calculation, the CUS model is adopted as the reference velocity model. Instead of specifying the focal mechanism for each event, here, the focal mechanisms are classified into three types (strike-slip for rake∼0°, dip-slip for rake∼90° and oblique-slip for rake∼45°) and the regression relation between the focal depth and the optimal period is calculated for the three cases.

Stable amplitude spectra for 18 events are obtained with the do_mft analysis, while for the other three events, multiple reverberations in sediments are too strong to identify short-period Rayleigh waves. Then, we determine the focal depth of the 18 events following the procedures in Fig. 9(b). For most events, the difference between the resolved depth (Depth_OP) and the depth in the SLU catalogue (Depth_SLU) is within 2 km (Fig. 14b and Table S5, Supporting Information). Excluding the two outliers (event nos 8 and 15), the average difference is −0.5 km and the STD is 1.6 km. For event no. 8, the moment magnitude is 3.6, and there is only 1 local station (distance of 95 km). The large difference between the Depth_OP (11.1 km) and the Depth_SLU (5 km) may be due to the low SNR of the waveform. There are two events (event nos 2 and 17) close to event no. 8, and the depths are 2.7 and 19.6 km, respectively, making it difficult to distinguish whether event no. 8 is a shallow earthquake or not. For event no. 15, the moment magnitude is 3.8, and there is only 1 local station (distance of 70 km). Although the Depth_OP (7.1 km) is much shallower than the Depth_SLU (18 km), it is consistent with the USGS regional moment tensor solution (7 km) and maybe more data are required to distinguish the accurate depth. Overall, even though only the simplified focal mechanism is used, the focal depth could be determined stably for small earthquakes.

Forward tests and applications to many earthquakes suggest that this method works well for shallow events (≤20 km), and here we investigate whether this method could be applied to deeper events (20–40 km). The method is based on the nearly linear relation between the optimal period and focal depth. Thus, we calculate the theoretical optimal period for deeper events based on eigenfunctions. The same four-layer velocity model (Table S1, Supporting Information) is used, and we take the vertical strike-slip faulting (0°/90°/0°) with azimuth of 45° as an example to calculate the surface wave amplitude spectra. The result shows that the optimal period for depth 20–40 km also increases with depth near-linearly (Fig. S11, Supporting Information), but the slope is significantly larger than the shallower cases, which may be due to the faster wave velocity in the lower crust. This linear relationship remains for depth of 20–40 km, suggesting the possibility to apply to deeper events. However, for deeper events, the minimum epicentre distance for observable Rayleigh wave is larger. For example, if the source depth is 30 km, the minimum epicentral distance is about 150 km, making it less feasible to apply to small earthquakes.

There may be complexity in the Rayleigh wave amplitude spectra at certain distances or regions. For example, the synthetic test based on eigenfunctions (Fig. 1a) shows that, for an extremely shallow event (∼ 1 km), the optimal period is ∼0.7 s and there is another local maximum at ∼5.6 s. If the optimal period is shifted towards a longer period because of high attenuation, there may be the possibility to take the local maximum at ∼5.6 as the optimal period in the automatical measurement, especially at farther distances (e.g. ∼70 km). The consistency between optimal periods measured at different distances provides a criterion to distinguish which is the unbiased optimal period, and a careful visual inspection at close stations may be helpful.

In this method, 1-D velocity model is usually adopted as the reference velocity model and lateral variation has not been considered. For regions with strong 3-D heterogeneities, such as the Los Angeles region, both the phases and the amplitudes of the observed Rayleigh waves may differ significantly from the predictions by the 1-D velocity model (Wang & Zhan 2020), and the scattering, focusing or defocusing may lead to unneglected errors or possible failure of this method, especially for very shallow earthquakes. We take the 2020 south California earthquake as an example and explore the impact of 3-D velocity model. We use the SPECFEM3D algorithm (Tromp et al. 2008) to perform simulations and the SCEC Community Velocity Model version 4 (CVM-4, Kohler et al. 2003) is adopted as the reference seismic velocity model. The realistic surface topography is interpolated from SRTM15 (Tozer et al. 2019). The modelled region ranges from 116.5° to 119.5° W and 33° to 35.5° N, which are divided into 288 × 288 grids, and the minimum observation period of the synthesized waveforms is ∼ 1 s. The focal depth (16 km) and other source parameters are from the USGS catalogue and regional moment tensor solution (Table 1). Then, the 1-D model is derived as an average from the CVM-4 model, based on which we calculate the synthetics as comparison. Then, we measure the optimal period of the 3-D and 1-D synthetics and calculate the difference at each station. The differential optimal period ranges from −0.1 to 2.0 s with the average value of 0.7 s (∼8 per cent of the average optimal period). We investigate the relationship between the differential optimal period and the tomography and topography (Fig. S12, Supporting Information), and have not observed obvious correlations between the differential optimal period and tomography (i.e. Vs at the depth of source or ground surface), but high elevation roughly corresponds to larger optimal period, which is consistent with the observation of real data (Fig. 12 ). Detailed analysis on the effect of small-scale velocity anomaly and topography would be explored in future works.

6 CONCLUSION

In this paper, we find that the optimal period of short-period Rayleigh waves at local distances is almost linearly correlated to the focal depth. Thereafter, we propose a new method to resolve the focal depth of small to moderate earthquakes automatically with the optimal period of Rayleigh wave amplitude spectra. Synthetic tests suggest that the method is robust against different source parameters and crustal structure. Applications to tens of Mw3.5–5.0 earthquakes in different tectonic settings suggest the method could be used to obtain reliable focal depth. In the central and eastern U.S., the average difference is about 0.5 km between the depths determined by the new method and the SLU catalogue.

In this method, the focal depth is obtained straightforwardly based on the observed optimal period and the near-linear relationship and does not require complex waveform data inversion, without need for precise focal mechanism. This advantage makes the method an useful tool to automatically and rapidly determine the focal depths of small to moderate earthquakes. Focal depth can be robustly determined even with one or a few seismic stations, implying applicability in sparse networks or historical events. Although the method has the advantage of being straightforward and automatic, combination of more observables (such as spectral null of Rayleigh waves, or P-wave polarization, etc) would lead to even more robust determination of focal depths, which can be explored in future works.

SUPPORTING INFORMATION

Figure S1. This shows the variation of eigenfunctions against focal depths.

Figure S2. This shows the non-normalized Rayleigh amplitude spectra.

Figure S3. This shows the senstivity test on the uncertainty of the focal mechanism.

Figure S4. This shows the senstivity test on a basin-range crustal model (WUS) and a stable shield model (CUS).

Figure S5. This shows the sensitive test on the non-DC component.

Figure S6. This shows the forward test on the applicability of the regression relation to events with different faulting mechanisms.

Figure S7. This shows more information and a histogram for the application to the 2014 Oklahoma event.

Figure S8. This shows the forward test on how the averaging scheme works.

Figure S9. This shows the application to the 2019 M_w 4.0–4.3 Weiyuan sequence

Figure S10. This sshows the application to the 2022 M_w3.7 Weifang earthquake.

Figure S11. This shows the relationship between the focal depth and the optimal period for deeper events (20–40 km).

Figure S12. This shows the investigation on the impact of 3-D velocity model.

Table S1. Velocity model used in the synthetic test. The four-layer model is based on the AK135 global velocity model (Kennett et al. 1995) and the CUS model (Herrmann et al. 2011).

Table S2. The CUS Velocity model (Herrmann et al. 2011).

Table S3. Source information for 2019 M_w 4.0-4.3 Weiyuan earthquake sequence. The fault solutions are from Yang et al. (2020).

Table S4. Source information for the 2022 M_w  3.7 Weifang earthquake. The fault solutions are from the waveform inversion in this study.

Table S5. Source information of the 21 M3.5+ events in the central and eastern U.S. (2014–2021). The point source solutions are from the St. Louis University (SLU) Earthquake Center.

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

This work made use of the Generic Mapping Tools (GMT, Wessel et al. 2019) and Seismic Analysis Code (SAC). Waveform data used in this study were downloaded from the Data Management Center of the Incorporated Research Institutions for Seismology (http://ds.iris.edu/ds/nodes/dmc/, last accessed 2022 May) and Data Management Centre of China National Seismic Network at Institute of Geophysics (SEISDMC, doi: 10.11998/SeisDmc/SN). This study is supported by The Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (2019QZKK0901), National Natural Science Foundation of China (grant nos 41804039 and 42030301), Guangdong Province Introduced Innovative R&D Team of Geological Processes and Natural Disasters around the South China Sea (2016ZT06N331) and Shanghai Sheshan National Geophysical Observatory (no. SSOP202106). We thank Yi Wang and Jun Xie for helpful comments. We thank editor Gabi Laske and two anonymous reviewers for constructive comments that helped to improve the paper. Author contribution statement: Xiaohui He: Methodology, Conceptualization, Writing original draft. Peizhen Zhang: Review & Editing. Sidao Ni: Methodology, Review & Editing. Risheng Chu: Review & Editing. Wenbo Wu: Methodology, Review & Editing. Kaiyue Zheng: Software.

DATA AVAILABILITY

The data underlying this paper will be shared on reasonable request to the corresponding author (e-mail [email protected]).

References