1 Introduction

Focal depth of seismic sources is one of the most important parameters in seismological studies. Accurate determination of depth is crucial in the study of seismo-genesis and seismic hazard estimation. Due to distributed nature of earthquake rupture slip, focal depth may be defined as the depth of rupture initiation point (hypocentral depth) or as the averaged depth weighted by slip (centroid depth). Hypocentral depth is usually determined with arrival times of P or S waves, while centroid depth is resolved from waveform information. For the purpose of rapid seismic hazard assessment (for example, theoretic prediction of shakemap), centroid depth is the essential parameter. But, the routine earthquake catalogs by NEIC, USGS CMT or Global CMT mostly involve very long period waves; thus, the centroid depth uncertainty is large, especially for shallow earthquakes, and further studies are needed to improve the depth accuracy (Engdahl et al. 1998).

On April 20, 2013, an Mw6.6 earthquake occurred near the border between Tibetan plateau and Sichuan Basin (the event is referred to as Lushan earthquake hereafter) (Fig. 1), and it caused substantial fatality and economic loss. Within 1–2 h of the earthquake, theoretic shakemap was constructed by CEA and helped emergent rescue. But, construction of shakemap requires reliable earthquake source parameters including centroid depth. After the Lushan earthquake, various groups provided different estimates of centroid depths, varying from 12 to 22 km (Liu et al. 2013; Shan et al. 2013; Xu et al. 2013). The difference in depth estimate is substantial, because a centroid depth of 22 km would predict much weaker ground damage as compared to depth of 12 km. Thus, it is worthwhile to resolve the accurate centroid depth and investigate factors that might influence centroid depth inversions. Although the difference in depth inversion may be due to different dataset used (local waveform data vs. teleseismic data, strong motion records vs. broadband velocity seismic records, etc.), here, we examine the issue of centroid depth inversion with only teleseimic body waves which are very sensitive to centroid depths and less suffer from crustal heterogeneity in the source region (Kikuchi and Kanamori 1982).

Fig. 1
figure 1

The Lushan earthquake (star) and global seismographic network stations (triangles). Stations from network AK and AU are also used in this study, but not plotted

Full size image

Teleseismic cut and paste (tCAP) method utilizes teleseismic (epicenter-station distance from 30° to 90°) body wave to invert source parameters (seismic moment, fault plane solution, and centroid depth) by assuming that earthquake source is a double couple (Chen et al. 2012). It is believed to be reliable for depth estimation for this method exploits depth phases (pP, sP) (Zhan et al. 2012). However, the dependency of the centroid depth inversion has not been systematically investigated regarding factors such as station azimuth gap, filter corner frequency, and source time duration. Here, we check the probable factors which might influence derived centroid depth and investigate how serious these factors can affect the inversion.

2 Data and method

2.1 The tCAP method

The original cut and paste (CAP) method was proposed to study earthquake source parameters for local waveforms (typically with epicentral distances less than 500 km) by Zhao and Helmberger (1994), Zhu and Helmberger (1996) in order to avoid the dominance of surface wave in the full waveform inversion, as well as to account for inaccurate velocity models via shifting segments of Pnl and surface waves separately. In the local CAP method, centroid depth, fault plane solutions (strike/dip/rake) are inverted with grid-searching algorithm, thus capable of finding the global minimum. Therefore, the local CAP method has been broadly applied in source parameter inversions for earthquakes in US, China, Iran, and other regions (Huang et al. 2009; Wei et al. 2009; Tan et al. 2010). The similar idea in local CAP method was applied for teleseismic body inversion to account for usually much stronger SH waves as compared to P waves, and the method is called as tCAP method (Wei 2009). Just as the local CAP, source parameters such as centroid depth and fault plane parameters are also inverted with grid search. The tCAP method has been demonstrated to be effective for determining centroid depth of moderate earthquakes (for earthquakes around M6) due to sensitive constraints from depth phases of pP, sP for P waves, and sS for SH waves (Wei 2009; Zhan et al. 2012).

2.2 Data

The Lushan earthquake is well recorded by global stations for both teleseismic P and SH waves (Fig. 1). We obtain three component broadband waveform data from the fast archive recovery method products of IRIS WILBER II system. The original data in seed format is converted to SAC format. Then, we remove linear trend and instrument response and rotate the horizontal components to the great-arc path. We pick the P/SH arrival time with theoretic time table from IASP91 model. We choose the data with signal-to-noise rate above 2, and noisy waveform data were not included for inversion. Time window of P and SH waves are set to be 45 s, long enough to contain complete P and SH wave trains, because the earthquake rupture lasted less than 30 s (Wang et al. 2013). We choose 0.01–0.16 Hz band-pass filter for P-wave inversion while 0.01–0.08 Hz for SH wave.

We compute the teleseismic Green’s Functions with the combination of plane wave propagation matrix method in the source region and ray geometry spreading in the mantle (Kikuchi and Kanamori 1986). A 1D velocity model from Tian et al. (2013) is used for the source region. Green’s Functions are calculated for depth range of 2–25 km, with interval of 1 km. Rupture duration is set to be 6 s, which will be discussed in the following sections. Then, for each depth, we do the grid search on source parameters (seismic moment, fault plane solution) for the best fitting between observed and synthetic waveforms computed from the Green Functions. The optimal centroid depth is claimed when the smallest waveform mismatch is achieved (Fig. 2a). In order to determine the optimal depth more accurately, we fit the mismatch-depth curve with quadratic fitting (Fig. 2b). The lowest misfit corresponds with the depth of 12.5 km. Fault plane parameters (strike/dip/slip) are 205°/47°/87° (Plane 1, or FM1) and 29°/43°/93° (Plane 2, or FM2). Moment magnitude (Mw) is 6.6. The fault plane solution and Mw results obtained in this study are consistent with catalogs of NEIC, Global CMT and with previous studies, arguing for effectiveness of tCAP method (Zeng et al. 2013; Xie et al. 2013; Han 2013). Though the inverted centroid depth of 12.5 km is very different from some previous results, it agrees with the results by Zeng et al. (2013) and (Han 2013). We defined the depth of 12.5 km as the optimal centroid depth of Lushan earthquake and defined uncertainty of centroid depth as the difference between this optimal depth and the inverted centroid depth in the following section.

Fig. 2
figure 2

tCAP modeling for Lushan earthquake. a Waveform fits. Blacktraces show observed seismogram, and gray lines are synthetic seismograms. Numbers to the left of the seismograms are time shifts (upper) and cross-correlation coefficient in percentage (lower). Positive time shifts indicate that synthetic waveforms have been delayed. b Waveform mismatch versus centroid depth

Full size image

2.3 Analysis of inversion parameters for centroid depth

It’s widely believed that several factors could influence the result of seismic waveform inversion: azimuth gap of station distribution, distance weighting, and the choice of band filtering. In order to investigate how these factors affect the inverted centroid depth, we perform tCAP inversion in the following sections.

Firstly, we study the uncertainties of centroid depth due to azimuthal distribution of seismic stations. We present the azimuthal gap effects in two ways. In the first way, we invert centroid depth for all stations only in a limited range of azimuth, and this is equivalent to maximum azimuthal gap. We divide azimuth range of 0°–360° into n sectors and invert centroid depth with seismic stations only in one of the sectors. We study the case of n from 2 to 6 (Fig. 3). That is also to say, maximum azimuth gap is from 180° (n = 2), to 240° (n = 3), to 270° (n = 4), to 288° (n = 5), and to 300° (n = 6). The centroid depth varies considerably (±3 km) even for n = 2. And, the uncertainty of centroid depth is larger for larger n, which is up to about 6 km when n equals 6. This implies that even with enough number of stations, azimuth coverage still substantially affects the result of centroid depth inversion.

Fig. 3
figure 3

Depth inversion results versus station azimuthal coverage (in terms of maximum azimuthal gap). The entire azimuth of 0°–360° is divided into n sectors (n = 2, 3, 4, 5, 6; from bottom to top, left to right). The outer-rings of the rose diagrams indicate the depth inverted with all stations in that sector

Full size image

We also test the azimuthal effects of minimum azimuthal gap. In this case, the entire station azimuth coverage is from 0° to 360°, but the azimuthal difference of neighboring stations has to be larger than specified value (the minimum azimuthal gap). Unlike the case of maximum azimuthal gap where many seismic stations are available for inversion, only limited number of seismic stations is available for the case of minimum azimuthal gap. We study minimum azimuthal gap of 15°, 30°, 45°, 60°, and 75° (Fig. 4). Intuitively, the uncertainty of centroid depth is expected to be higher with larger minimum azimuthal gap. However, since the distribution in this case is much more average, the derived depth should not vary too far. Result of the centroid depth inversion in Fig. 4 indeed demonstrates the expectation. Note that number of seismic stations involved in inversion drops significantly as the minimum azimuth gap goes up.

Fig. 4
figure 4

Depth inversion results versus minimum azimuthal gap. The azimuthal gap increases from 15° to 75°

Full size image

Then, we proceed to investigate the effects of exclusion of SH waves in centroid depth inversion. Previous studies involving tCAP method use both P and SH body waves (Xie et al. 2013; Wei 2009). For moderate events (M6 or weaker), SH waveform is typically noisier than P waves, and that is probably why some studies only use P waves for tCAP inversions (Zeng et al. 2013). Here, we exclude SH waves in the inversion process and find that the centroid depth for the case P wave only is 12.0 km, only 0.5 km different from the case of using both P and SH waves. Thus, SH waveforms do not seem to affect centroid depth inversion substantially. In the meantime, the strike is about 20° different from the preferred mechanism, which is mainly because teleseismic P-wave radiation for a thrust mechanism on 45° dipping fault has the same polarity and thus not sensitive to the strike (Fig. 5). Therefore, it is important to include SH waves to better resolve the mechanism, although we mainly focus on the depth resolution in this study.

Fig. 5
figure 5

tCAP modeling for Lushan earthquake

Full size image

Then, we explore the effects of source time duration (or alternatively, the rupture duration plus rise time) on centroid depth inversions. Usually, shorter source time duration is helpful for more accurate centroid depth inversion, because the depth phases (sP or pP, or sS) can be sufficiently separated from direct phase of P and SH. In some previous studies (Zeng et al. 2013; Wei 2009), the source duration is given as a priori parameter based on the scaling law with seismic moment magnitude (Dreger and Helmberger 1993; Ichinose et al. 1998). However, in fact, source time duration can vary remarkably from the predicted value. We test the effects of source time duration for two scenarios. In the first case, the centroid depth is fixed at 12.5 km, and source time duration is varied from 1 to 10 s with interval of 1 s so as to compute waveform mismatch. The misfit curve in Fig. 6 shows the best source duration of 6 s. In the second case, we invert centroid depth for each duration in the range of 1–10 s with Green’s Function for each depth grid. From Fig. 7, it is observed that different depth and source duration can produce the same waveform misfit which implies the trade-off between these two parameters. Thus, the centroid depth inversion substantially depends on source parameters. The inverted depth varies from 12.5 to 22 km, large enough to explain the range of depths from various studies (Global CMT; USGS NEIC CMT; Xie et al. 2013; Zeng et al. 2013; Han 2013) Fig. 8.

Fig. 6
figure 6

Waveform misfit versus source duration. Number above the beach-ball is the corresponding Mw

Full size image
Fig. 7
figure 7

Inverted depth versus different source duration. The dash line indicates the depth of 12.5 km (from Fig. 2). The black curve indicates the best depth obtained for different source durations. The gray line shows the waveform mismatch at the best depth for different source durations

Full size image
Fig. 8
figure 8

Inverted depth versus corner frequency of low-pass filter. Each star represents the depth inverted for the corner frequency

Full size image

We then explore the effects on centroid depth inversion from different band-pass filters with fixed lower cut off frequency in the inversion process. Usually, it is assumed that long enough period low-pass filtering is good enough for waveform inversion (the filter period has to be longer than the source time duration of the earthquake, to make sure that the point source assumption is valid). But, when the filter period is too long, the waveform information regarding depth phases may be filtered out, thus losing depth sensitivity. Our tests on different low-pass filter periods show that the results of centroid depth inversion are stable as long as the corner frequency of the low-pass filter is not too low. The centroid depth varies about 1–2 km when the corner frequency goes up, which is relatively smaller compared with the uncertainty due to azimuthal gap or source durations. However, the centroid depth in the waveform inversion with 0.06 Hz (or lower frequency) filter corner frequency is around 17 km, implying that choosing excessively long period will lead to biased centroid depth.

Velocity model also plays an important role in centroid depth inversion as the travel time difference between direct phase and depth phase can be fitted with different depth and velocity combinations. Here, we tested another 1D velocity model extracted from CRUST2.0 at the epicenter (Bassin 2000), which turns out to be an extremely low near surface velocity model (Fig. 9a). The depth obtained by the CRUST2.0 model is shown in Fig. 9b in which the depth estimate drops to 11.3 km. This value is 1.2 km different from the optimal centroid depth. Considering the model applied here which is an extreme low-velocity basin structure, this is acceptable. Note that the time shift allowed in the CAP algorithm can partly reduce the influence introduced by the velocity model uncertainty.

Fig. 9
figure 9

Depth inversion test under different velocity models. a two velocity models used in this study. Solid lines represent the model used in previous figures, while dash lines represent the sedimentary structure extracted from CRUST2.0. b Waveform mismatch versus centroid depth for the inversion using CRUST2.0 model

Full size image

3 Discussion and conclusions

We apply tCAP method for global broadband waveforms to invert centroid depth of 2013 Lushan earthquakes. We apply a series of tests on possible factors which could affect the inversion results of depth. We find that the key for reliable depth estimation from global teleseimic waveform is that the station azimuthal coverage should be as even as possible, and the number of stations should also be guaranteed (8 or more). This requirement of reliable centroid depth is higher for inverting fault mechanism parameters where only six stations are needed to get stable results (Chen et al. 2012).

The source duration choice should be rather careful, because the inverted depth is sensitive to source time duration. A feasible approach is to check the misfit-source duration curve for each depth grid. We also find that very long period band pass makes the depth inversion worse. Very long period band pass is widely used in waveform inversion for popular catalogs, such as GCMT and USGS. But, the depth information is lost during the low-frequency filter process; this could explain why most depth results in these catalogs need to be re-examined with shorter-period filtering.

In this paper, we discuss several influencing factors on the depth inversion, taking Lushan earthquake for example. However, we assume 1D velocity structure in this study and have not explored the error of centroid depth due to 3D crustal and mantle structure. To get more comprehensive analysis on uncertainties of centroid depths, future studies involving 3D waveform modeling are necessary.