Accurate Redetermination of the Focal Depth by Using the Time Intervals between the Inner Core Phases PKIKP and pPKIKP

Abstract

The hypocenter parameters of an earthquake may give us an insight into the Earth’s structure and tectonic processes. Among the hypocenter parameters, the focal depth is normally more difficult to estimate than the earthquake location (latitude and longitude). We propose to use the pPKIKP-PKIKP arrival time intervals for estimating the focal depth. We analyze the sensitivity of the pPKIKP-PKIKP arrival time interval to the earthquake depth. We measure the pPKIKP-PKIKP arrival time interval on seismograms (the vertical component), and invert the time interval data set using the simulated annealing inversion algorithm. We illustrate the inversion approach on two teleseismic earthquakes which have shallow and deep focal depths, and demonstrate that the approach is indeed appropriate to the shallow and deep event. We can obtain a reliable estimate on focal depth, even though the seismic station is sparse or in a remote part of the epicenter.

1. Introduction

Hypocenter parameters, such as latitude, longitude, depth, and origin time, are a crucial parameter in seismology, and the research work on the field of seismology is particular pertinent to hypocenter parameters [1,2,3,4,5]. Using the time interval between the arrival of the primary and secondary waves (S-P time interval) that is measured at each station, earthquake location, i.e., latitude and longitude, can be used to determine accurately and easily by analyzing the S-P time interval. However, accurate measurement of the focal depth of an earthquake is a typically more challenging and rewarding work, as the source depth is usually inadequately constrained [6,7]. Therefore, precise determination of source depth is normally more difficult than determination of location, unless denser seismometers are deployed in the vicinity of the earthquake epicenter. Furthermore, it is well known that precise depth determination for shallow teleseismic events is a very complicated subject [8], and the S-P arrival time interval becomes less and less sensitive to hypocenter parameters as the epicenter distance increases [7,9].

Seismologists have long recognized that high precision of focal depth remains highly problematic in seismology [10]. Craig [11] gave a general overview of the existing methods to estimate focal depth, which mainly contain waveform-based methods and traveltime-based methods. Heyburn et al. [12] estimated focal depth by measuring both teleseismic depth phases pP and sP and surface wave amplitude spectra. Florez and Prieto [10] proposed an array processing technique to calculate focal depth using pP-P times and double-difference algorithm, which needs two close events. Craig [11] presented a method for accurate depth determination for moderate magnitude earthquakes using teleseismic depth phases pP, sP, and sS data sets. Yuan et al. [7] developed Depth-Scanning Algorithm to obtain focal depth at local and regional distances, which was based on the traveltime residuals between the direct and depth phases (Pg, pPg, sPg, Sg, sSg, PmP, pPmP, sPmP, SmS, sSmS, etc.). Among all the methods mentioned above, the data used for determining focal depth are from the phases penetrating above the mantle-core boundary, and the methods can not be well constrained. However, the teleseismic depth phases, which almost travel across the whole earth, are well known to assess focal depth and can be a clear improvement on the depth determination accuracy [5,10,13,14]. Since the teleseismic depth phases are sensitive to hypocenter depth, they can be used to obtain more precise focal depth.

Firstly, we pre-process the raw seismograms from a seismic array deployed at Guangxi, China, and identify the arrival times of PKIKP and pPKIKP phases based on the peaks of waveforms and cross-correlation. Secondly, we invert the data set of pPKIKP-PKIKP time interval for the focal depth, using the simulated annealing (SA) method. Finally, we verify the accuracy and feasibility of the procedure using the teleseismic data from two events. These public domain data sets are from China Earthquake Administration (CEA), German Research Center for Geosciences (GFZ), International Seismological Center Engdahl-van der Hilst-Buland (ISC-EHB) Bulletin [15], and United States Geological Survey (USGS) catalogue.

Thus far, the time difference between the inner core phases pPKIKP-PKIKP has not used to determine earthquake depth, This may be because the arrivals of inner core phases are not clear-cut. On the other hand, due to a long travel through inner core with a nearly vertical incidence to the surface and losing most high frequency compositions, the inner core phases pPKIKP and PKIKP may be difficult to acquire from a teleseismic earthquake. However, if we choose an appropriate band filter for the teleseismic phases, pPKIKP and PKIKP should be identified on the seismograms more simply. Meanwhile, the phase pPKIKP is of a small incident angle, so that the phases pPKIKP almost take the same path as PKIKP between the source and station, which will prevent the effect from the lateral variation of velocity structure beneath the epicenter. Furthermore, with the increasing deployment of broadband seismometers worldwide, we can fetch more high-quality seismograms from a teleseismic earthquake, which will provide an abundant database to improve earthquake depth accuracy during the determination procedure.

2. Methods

2.1. PKIKP and pPKIKP Phases

Figure 1a demonstrates the seismic ray path of PKIKP and pPKIKP phases. Inner core phase PKIKP is a compressional refraction wave across downward the inner-core boundary (ICB), and pPKIKP is the phase which goes up firstly from the source and then is reflected at the earth’s surface downward to the same like PKIKP. To illustrate the sensitivity of the pPKIKP-PKIKP time interval and focal depth, it is assumed that the great circle from station to event is 150° (Figure 1a), and the hypocenter depth ranges from 50 to 600 km. The time difference pPKIKP-PKIKP is calculated using ak135 Earth model [16] and plotted as function of focal depth (Figure 1b). The slope of the function is about 0.22 s/km, thus the pPKIKP-PKIKP phase can offer a good constraint on focal depth. Compared to the theoretical pP-P time (Figure 1a in [10]), pPKIKP-PKIKP time differences are more sensitive to earthquake depth. In addition, the ray path in Figure 1 indicates that the pPKIKP-PKIKP time interval data is mainly affected by the velocity structure above the focal point.

2.2. Redetermination Algorithm for Focal Depth

To redetermine the focal depth, d, the arrival times of pPKIKP and PKIKP at the same station of seismic array will be the observation data sets. The focal depth is assumed at a location that best fits the observation and has the smallest misfit. The fit is resolved by comparing the theoretical and observed pPKIKP-PKIKP time interval.

Let us consider an earthquake which was recorded by seismic array and the observed pPKIKP-PKIKP time interval at ith station, referred to as $\Delta t_{obs}^i$

$$\Delta t_{obs}^i={(\mathrm{pPKIKP}-\mathrm{PKIKP})}_{\mathit{i}}$$

(1)

and $\Delta t_{cal}^i\left(d\right)$ represents the theoretical pPKIKP-PKIKP time interval calculated by ak135 model at a depth d for the ith station, thus the optimum focal depth d, $F\left(d\right)$ can be obtained by minimizing the time residuals between observed and calculated pPKIKP-PKIKP time intervals

$$F\left(d\right)=\sum _i{(\Delta t_{cal}^i\left(d\right)-\Delta t_{obs}^i)}^2$$

(2)

During the redetermination of focal depth, the earthquake epicenter, latitude, and longitude were kept the same as that in the ISC-EHB catalog, that is, the epicenter does not need to improve. Finally, the focal depth was redetermined by the global optimization method based on simulated annealing that minimizes the misfit Equation (2) [17,18].

2.3. Simulated Annealing

Since Kirkpatrick et al. [17] proposed the simulated annealing algorithm based on statistical physics, it has become an important tool for resolving nonlinear inversion issues. In this paper, we utilize the Metropolis algorithm to find the global optimum solution of Equation (2). Metropolis algorithm is a fixed temperature for the classical general scheme of the simulated annealing. In the inversion procedure, we use a temperature of 100 °C.

To numerically find the globally optimal solution of the objective function (2), we assume the model parameter m is the only variable to be inverted, $\mathrm{m}=d$. The model space m has a variation of plus or minus 30 km of the ISC-EHB depth for every iteration in the inversion. The model space must be greater than 0. Here are the instructions for how to do the simulated annealing procedure using the classical scheme.

• Basic variables are declared and set to their default values. In the inversion technique, the following parameters were set: (1) An initial model parameter $m_0$ was generated randomly, (2) the initial temperature $T=T_0$ is set to 100 °C, (3) The maximum number of iterations in the inversion procedure, (4) fitting error threshold $ϵ$.
• Estimate the model parameter ${\mathrm{m}}_k$ in the objective function (2), and the temperature $T_k$, for the kth iteration. The acceptance of the accepting model parameter ${\mathrm{m}}_{k+1}$ from the current ${\mathrm{m}}_k$ is given by the following transition probability:

  $$P(\Delta F)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\Delta F<0\hfill \\ \exp \left(\frac{\Delta F}{T}\right)\hfill & \mathrm{if}\Delta F\ge 0\hfill \end{array}\right.$$

  (3)

  where $\Delta F$ is the difference of the objective function (2), that is, $\Delta F=F\left({\mathrm{m}}_{k+1}\right)-F\left({\mathrm{m}}_k\right)$. The Boltzmann constant is not included in Equation (3). If $\Delta F<0$ then the transition probability is equal to 1, indicating that the system energy tries to go to the lower energy states. Thus, the new model parameter is accepted. However, if $\Delta F\ge 0$ then we should calculate the probability function $P(\Delta F)$ (The value of $P(\Delta F)$ takes a value from 0 to 1), and then a random number R between 0 and 1 is generated. If $P(\Delta F)>\mathrm{R}$ then the model parameter will be updated to ${\mathrm{m}}_{k+1}$, otherwise the ${\mathrm{m}}_{k+1}$ will be rejected.
• Evaluate the objective function $F\left({\mathrm{m}}_k\right)$. If $F\left({\mathrm{m}}_k\right)<ϵ$ then the inversion procedure will be terminated, otherwise, $k=k+1$, going to step 2.

3. Applications

There are two data sets that were used for the applications of the presented method. The first one is from a shallow crustal Earthquake that occurred on the 23 August 2014 Mw 6.2, Chile, and the second one is from a deep Earthquake that occurred on the 21 February 2017 Mw 6.5, Bolivia. The performance to redetermine the focal depth for these two data sets will illustrate that the method can redetermine the focal depth accurately. All data have the same sampling rate of 0.01 s. After data pre-processing for the raw data (mean and trend removing, filtering, picking etc.), the inner core phase PKIKP and pPKIKP were well identified for the performances of this method, and the redetermined focal depth will be compared to the results from four different agency, CEA, GFZ, ISC-EHB, and USGS.

3.1. The Mw 6.1 Chile Shallow Crustal Earthquake

A strong shallow crustal earthquake struck 23 km WNW of Hacienda La Calera, Chile (Figure 2) on 23 August 2014 at 22:32 UTC. According to the CEA, GFZ, ISC-EHB, and USGS earthquake catalogue, the event occurred at a depth range from 30 to 42 km (

Table 1. The result and comparison for the shallow event occurred on 08/23/2014, Mw 6.1, Chile.

Event ID	Latitude (°)	Longitude (°)	Depth (km)	Magnitude	Author	Misfit (s${}^2$ )
null	−32.710	−71.390	30.0	6.5	CEA	21.0
gfz2014qngj	−32.641	−71.284	32.7	6.4	GFZ	3.0
usb000s5rc	−32.695	−71.442	32.0	6.4	USGS	6.1
610572067	−32.688	−71.351	42.0	6.2	ISC-EHB	46.0
null	−32.688	−71.351	34.3	6.2	This method	0.1

) with less agreement. The seismic station is shown by a black triangle in Figure 2. The vertical component seismograms were band-pass filtered in the frequency band 0.1–0.6 Hz, normalized with respect to the maximum amplitude, and sorted by the distances. We have collected a set of 12 records with clear PKIKP and pPKIKP phases at distances ranging from 169.8° to 170.6° (Figure 2). The theoretic and observed traveltime of PKIKP are marked with red plus sign and filled circle, respectively, in Figure 2, and the blue indicates pPKIKP.

In the inversion procedure, we chose an initial depth between minus and plus 20 km with the depth of ISC-EHB at random. The depth estimation of this event by the pPKIKP-PKIKP time intervals was shown in Figure 2 and four sets of criteria by different bulletin are shown in Table 1. In this case, taking the focal depth in ISC-EHB as the initial value, the depth estimation by using pPKIKP-PKIKP time intervals is 34.3 km, and the residual misfit is about 0.1 s${}^2$, which is minimal in comparison to the four catalogs. The result also indicated that this method can estimate the focal depth for shallow events by using teleseismic pPKIKP and PKIKP phases.

3.2. The Mw 6.5 Bolivia Deep Earthquake

This earthquake measuring Mw 6.5 struck Southern Bolivia on 21 February 2017. According to the event catalogs of four differences agency bulletin, the focal depth is between 593 to 602 km (

Table 2. The result and comparison for the deep event occurred on 02/21/2017, Mw 6.1, Bolivia.

Event ID	Latitude (°)	Longitude (°)	Depth (km)	Magnitude	Author	Misfit (s${}^2$)
null	−19.170	−63.800	600.0	6.3	CEA	47.0
gfz2014qngj	−19.264	−63.931	592.9	6.5	GFZ	164.8
usb000s5rc	−19.281	−63.905	596.0	6.5	USGS	104.7
610572067	−19.296	−63.956	602.5	6.5	ISC-EHB	22.5
null	−19.296	−63.956	608.0	6.5	This method	0.8

), which is a typical example of deep Earthquake. After pre-processing and filtering by passband of 0.05–0.3 Hz, the teleseismic data (vertical component) from 18 stations with the distances ranging from 170.4° to 174.5° were selected to yield pPKIKP and PKIKP data sets (Figure 3). Perhaps due to a time synchronization delay during data acquisition, a larger time delay between theoretical and measured arrival time of pPKIKP and PKIKP phases presents at a distance of 172.5° (the magenta-colored waveform in Figure 3). However, such a case does not ruin the proposed method to produce an accurate measurement due to using pPKIKP-PKIKP time intervals.

Taking the initial depth as the focal depth in ISC-EHB with a deviation of ±30 km, the focal depth estimated for this event is 608.0 km, as shown in Table 2, and the residual misfit is only just 0.8 s${}^2$, which is the minimal among the results from different bulletins.

4. Conclusions

In this study, we present a straightforward approach to determine focal depth by utilizing the time interval between pPKIKP and PKIKP. Using teleseismic records acquired at South China with the event of the Chile shallow crustal earthquake and Bolivia deep earthquake, we confirmed that the presented method can accurately determine the depth of teleseismic earthquakes. The determining earthquake depth approach was extended to cover the inner core phase.

We showed that teleseismic pPKIKP-PKIKP time difference can be used to calculate the depth of an earthquake from shallow to deep. The pPKIKP-PKIKP time interval were more to the depth than to the depth phase, which means pPKIKP-PKIKP can constrain the depth of teleseismic earthquakes efficiency. In order to obtain the reliability of pPKIKP-PKIKP time interval, we suggest filtering the seismograms records with different bandwidth that can get a clear pick of those two phases.

The proposed new method, by using the observed data with the arrival time of pPKIKP minus PKIKP, can calculate earthquake depth precisely. By comparing the focal depths published in different catalogs and bulletins, the results from two typical applications indicate that the proposed method can estimate earthquake depth much more accurately, with a minimal misfit both for shallow and deep events.

The method we propose could be used to redetermine a precise focal depth by using the data from portable or permanent seismic array, no matter if the seismic station is really sparse or the station network is far away from the earthquake epicenter, and disregarding whether the time synchronization is correct when the data were acquired.