Investigation of Source Depths of the 2006 Pingtung Earthquake Sequence Using a Dense Array at Teleseismic Distances

We determined source depths of 6 major events of the 2006 Pingtung Earthquake sequence by comparing the observed array beams at teleseismic distances with synthetic seismograms for a range of depths. A simple but robust procedure to identify the depth of a seismic event has was developed and successfully demonstrated through its application to the Pingtung Earthquake sequence. This method is theoretically based on the fact that adopting the surface-reflected waveforms at teleseismic distances can eliminate the tradeoff relationship between the source depth and epicenter determination. Accordingly, we utilized high quality seismic records of the Yellowknife Array with our method and found that we did obtain reasonable results; the depths of 6 events previously determined by CWBSN were overestimated except a second large event. This event was a complex source rupture, which occurred on December 26, 12:34 UTC (Mw = 6.9), based upon a detailed examination of near source strong motion seismograms. Combining analyses from near source and teleseismic observations, we suggest that the shallow fault plane rupture progressed deeper during the event. Investigations based on teleseismic and local observations identified different portions of the rupture fault, wherein the CWBSN data constrained the initial onset and the Yellowknife Array resolved the centroid of the rupture.


INTRODUCTION
On 26 December 2006, two closely timed earthquakes with magnitudes of (M w ) 7.1 and 6.9 respectively occurred at intermediate depths offshore of southwestern Taiwan.The first earthquake started at 12:26 UTC (20:26 LT).According to the Central Weather Bureau (CWB) catalog, the epicenter was located at 21.67°N and 121.56°E, off the southwest coast of Taiwan, with a focal depth 44 km.It was followed by another earthquake eight minutes later at 21.97°N and 121.42°E (about 36 km to the north-northwest of the first shock), which had a focal depth of 50 km.Both events, known as the 2006 Pingtung offshore earthquake doublet, occurred within a zone of transition along the north-south boundary between the Eurasian plate and the Philippine Sea plate.The normal-faulting focal-mechanism of the first shock suggested that this shock occurred as the result of intraplate stresses within the subducting Eurasian plate.Historically, this region is low in seismicity.The twin earthquakes were the biggest events which have occurred to the west of Hengchuan Peninsula in the past 40 years.A total of 1774 aftershocks occurred within one month after the Pingtung earthquake as determined by the Central Weather Bureau Seismic Network (CWBSN), a short period seismic network composed of 75 stations entirely covering Taiwan and its off-shore islands to monitor regional seismicity (Fig. 1).However, almost all of stations of CWBSN are located on the islands and can not provide a ideal coverage for off-shore earthquakes such as the Pingtung earthquake sequence.This results in strong tradeoff estimations between epicenters and source depths of those events when only using observations from CWBSN.Specifically, source depths and the tectonic origin of the Pingtung earthquake sequence are still a matter of debate.An accurate determination of source depths of the Pingtung earthquake sequence is important for improving the accuracy of aftershock patterns as well as providing much needed information about the fault plane orientation.
By introducing surface-reflected phase arrivals at teleseismic distances, clear arrivals on the seismic record are theoretically represented for different paths and have been employed to identify source depth of earthquakes (Engdahl 2006).Theoretically, those surface-reflection P phases have similar ray paths as a first-arrival P waves; adopting teleseismic records could eliminate the tradeoff relationship between source depth and epicenter determination.Taking advantage of array observations at teleseismic distances, the stacked array waveform provides high resolution seismic phases for source depth estimation.In this study, the best depth solutions of 6 major events of the Pingtung Earthquake sequence are accurately obtained using array observations at an epicentral distance near 85°.In sum, considering the arrival times of surface-reflected phases for locating source depths, it will help us easily control the range of depth and improve accuracy.Furthermore, a reliable earthquake depth estimation can provide valuable data for improving our understanding of Taiwan orogeny.

DATA
After searching for available global seismic data, we used data from the Yellowknife array (YKA), a small-aperture array consisting of 18 short-period, vertical component instruments with a dominant response frequency of approximately 1 Hz, plus four broadband sites with special sensors which detect a wide range of seismic wave frequencies.This seismic array located near Yellowknife in Canada's Northwest Territories was installed in 1962.This seismological array is located nearly 85 degrees away from the Pingtung earthquake sequence and seismic P waves and its surface-reflection phases from this range can avoid interference by other complex reflections near the source region.Furthermore, the location of the Yellowknife far away from coastlines and sources of man-made noise guarantees Fig. 1.Distribution of aftershocks of the Pingtung earthquake sequence reported by CWBSN.Large dots show the locations of major events analyzed in this study.Each number indicates an event number in Table 1.
very low noise conditions and offers a good quality record.The array aperture is 20 km and the instruments are deployed along two perpendicular lines (in N-S and E-W directions) with an interstation spacing of 2.5 km (Fig. 2).The YKA was designed specifically to detect high-frequency P waves and is therefore especially well-suited to study the short-period seismic wavefield at regional and teleseismic distances (Weichert and Whitham 1969).The site of the YKA is located in archaean granite-greenstone block located in the northwestern Canadian shield.The structure of the crust beneath the array is highly homogeneous and horizontally layered (Corbishley 1970), and the velocity structure is well documented by a number of local studies.The layering beneath the array is horizontal so that little or no variations in azimuthal travel time exist (Hwang and Clayton 1991).YKA can be considered as an ideal dense and small array suitable for identification of earthquake source depths and rupture properties of an earthquake sequence at teleseismic distances.
For this study, we requested 6 large sequential events included the first event which occurred on 26 December 2006, 12:26 UTC.The epicenters, depths, magnitude (M L ), average source to array distances, and average azimuth are listed in Table 1.The epicentral distances of all events are larger than 85°and the azimuthal differences less than 0.1°.Though the epicenters of the selected 6 events are distributed in a small area as shown in Fig. 1, the variations are insignificant in comparison with the epicentral distance and azimuth.We can treat the paths of these events as the same and therefore we expect to observe similar phases for future analysis.Moreover, we can discount pwP at the water depths less than 1.5 km because it is nearly impossible to separate the pP and pwP arrivals theoretically for most records (Engdahl et al. 1998).

DATA ANALYSIS AND RESULTS
Rays traveling upward from the hypocenter to the free surface are called surface reflection phases or depth phases (pP, pwP, sP); among these, the pP phase is the one reflecting and traveling to the recording station as a compressional wave (Fig. 3).We can see the incident ray paths as nearly perpendicular to the source and station while the source to station distance is greater than 80° (Pho and Behe 1972), so the downward segment of pP is almost the same P (Fig. 3, right panel).Moreover, at distant stations, the differential time of pP to P phase changes slowly with distance but rapidly with depth and implies that the time difference of P and pP mainly result from the path difference between source and surface.In other words, the time difference between P and a depth phase is highly sensitive to the earthquake depth and unaffected by errors in estimates of origin time.Furthermore, referring to global travel time tables (Jeffreys and Bullen 1940;Kennett 1991;Kennett et al. 1991), at one station located at a distance greater than 80°, the arrival of PcP wave is very close to P and its signal even mixes with the P wave.This wave train can be separated from later pP plus pPcP wave trains.No other obvious arrivals between PcP and pP were expected.In considering the situation of point source, since the pP phase propagates a little bit longer path than P and dissipates more energy, so that the amplitude of pP phase is normally smaller than the direct P wave.Hence, we can pick the arrival of the secondary big amplitude as pP and recognize it easily in an observed seismogram from a station whose epicentral distance is greater than 80°.
To employ differential travel time information of P and pP to identify source depths, the primary step should clearly and correctly recognize the depth phase.In our study, we proposed to use slant stack procedures (Vidale and Benz 1992) to reduce the effects of random noise in the data and to enhance the depth phase signal.In fact, a stacked array beam did effectively reduce the noise and yield a better representation of the signal of interest than the individual seismograms.Employing this technique, we can clearly observe the pP in an array beam easier than in individual seismograms.Figure 4 shows an example of stacked P and pP arrivals from the array beam and the individual waveforms of the 4 th event listed in Table 1.The stacking process helps us to pick phase arrival time more accurately.For a deep earthquake, the depth phase is much later than the direct P wave, so it can be easily separated from the P phase.In contrast, for a shallow earthquake, surface reflections arrive shortly behind the direct arrival and may be interfered with by scattered P wave energy from both the source and receiver sides, so as to be hardly identified.
Actually, depth phases such as pP, sS from a shallow earthquake always arrive just behind the P phase and interfere with the P phase to produce a complex pulses.However, the shape of the interference packet changes with depth.This information supplies a means of estimating the source depth.Comparing observations with suitable synthetic seismograms (e.g., Langston and Helmberger 1975;Helmberger and Burdick 1979) can help us identify source depths accurately.
We adopted the Frequency-Wavenumber (FK) synthetic seismogram package of Zhu and Rivera (2002) to calculate synthetic waveforms.This is a forward modeling approach to simulate waveform using a 1D earth layer model, epicentral distance, source depth, azimuth and focal mechanism.It is well known that the travel time is highly correlated to the earth's structure, so the 1D velocity model we used here was consistent with the CWB Taiwan (Chen 1995) and IASPE91 1D models (Kennett 1991) for depths above 200 km and below 200 km, respectively.In addition, the focal mechanism will theoretically make the simulated waveform differ in relative amplitudes and polarities of P and pP waves (Murphy and Barker 2006) but it will not change the relative arrivals within those phases.Therefore, we selected all mechanisms reported by the global CMT catalogs of Harvard (GCMT), USGS, the regional CMT catalog of the Broadband Array in Taiwan for Seismology (BATS), and National Research Institute for Earth Science and Disaster Prevention (NIED) to compute synthetic seismograms, then selected the proper one which fit the observed records well.Herein, epicenters from the CWB catalog were employed.A different selection for the epicenter will slightly alter the absolute phase arrival but would not change the time delay between pP and direct P phase.
We calculated synthetic seismograms from one epicenter distance with different source depths following a constant depth interval to obtain a depth seismic profile (e.g., Fig. 5a).We can then move the beam trace on this seismic profile to find a suitable simulated waveform which fits the P and pP time difference well.Figure 5a shows a comparison of observed and synthetic waveforms for different depths for the 4 th event listed in Table 1.This shows that the beam trace was best fit by the synthetic waveform with a source depth of 34 km.The source depth search can be quantitatively represented by a cross-correlation analysis.The coefficients of the cross-correlation calculation within an array beam and synthetic seismograms provide a coherence map that presents the consistency of the time difference of P and pP phases within observed and calculated seismograms.The maximum value on the map means the minimum residual between the pP-P time differences between observed and synthetic records which corresponds to the depth value signify the optimum depth.It can help us to convert the observed pP-P time differences to the best fit source depths.Finally, we can infer the source depth for 4 th event as 34 km from the results of comparison of the stacked array beam and depth seismic profile and the cross-correlation calculation (Figs.5a, b).
The depths of those 6 events determined by CWBSN are  located between 40 to 50 km.Using our method, the depths are modified to fall between 32 to 56 km (Fig. 6), thus our results are slightly shallower than the CWBSN reported except for the 2 nd event listed in Table 2 (named Ev2 in this study).We determined the source depth of this event as 56 versus 50 km as determined by CWBSN.Notice that this event was the biggest event in the Pingtung earthquake se-quence reported by CWBSN and also is the only case where our new depth estimate is deeper than the CWBSN estimate.
To verify the complex source rupture of this event, several near source strong motion records were examined.Figure 7 shows 4 velocity seismograms from the nearest sites of Ev2.All of them clearly present direct P arrivals (marked by T1 in each seismogram).The predicted S ar-  rivals (marked by T2 in each seismogram) can be clearly recognized in the acceleration seismograms which arrive approximately 6 to 7 sec later than the P phases for all stations.However, in each seismogram, a complex vibration pattern follows the theoretical S arrival and implies that this event should not be taken as a simple point source.We found that large horizontal component amplitudes (marked by a red dashed line on each three component seismograms) appeared with a time delay of 8 -10 sec after the theoretical initial S arrival.This large pulse can be considered as a big rupture 8 -10 sec after initial rupture; however, we have no spatial resolution concerning this asperity based upon a limited number of near source observations.
Using the teleseismic waveform records shown in Fig. 8a, a complex waveform is seen in the stack beam (top trace of Fig. 8a).This is different from a typical aftershock waveform as shown in Fig. 4. We identified two emerging wave slots with a similar oscillation pattern and duration following the initial arrival; however, due to the synthetic program just simulated simple pulse waveform, emerge onset of this event restricts us to apply this method to determine source depth using depth phases as shown in Fig. 5. Therefore, we applied a low-pass filter with a signal less than 0.5 Hz on the array beam trace.After filtering, two clear wave trains were identified and it showed that this simple low frequency wave train arrived nearly 10 sec after its initial arrival (second trace of Fig. 8a).We considered that this low frequency wave train recorded seismic energy radiated from a large asperity.This wave train can correspond to that observed in a near source as shown in Fig. 7. Actually, two low frequency wave slots were clearly observed in this filtered seismogram.We consider that the first one is P and the second is pP, and both originated from the large asperity.The filtered trace was taken as the reference in our cross-correlation calculation.According to our processes, the estimated depth was determined by the time difference of this wave train and the estimated source depth should indicate the location of the large asperity and not the initial onset of earthquake fault.

DISCUSSION AND CONCLUSIONS
Traditionally, the arrival time of a predicted phase of a trial hypocenter and picked time of the observing seismogram are used for calculating the position of a seismic event.The optimal solution is obtained under the minimum phase arrival time residuals occurred (observed minus calculated).However, this approach is strongly correlated to the hypocenter and origin time.Especially for an event occurring outside of a seismic array, there is usually a strong tradeoff correlation between epicentral distance and origin time.Although we can remove the effect of the earthquake origin time by the use of differential times between phases (for example, S and P arrival) at the same station, there still exists a tradeoff relationship between source depth and distance range that leads us difficult to release the drawbacks.In suppressing the uncertainty, many studies (Engdahl and Fig. 6.Comparison of the source depths of analyzed events determined by the CWBSN (transverse axis) and by this study (vertical axis).The number represents the event number in Table 2. Gunst 1966;Engdahl et al. 1998;Wang and Zhao 2005;Warren and Shearer 2005;Engdahl 2006) have shown that using later arriving phase data can provide greater constraints on hypocenter parameters.When depth phase arrival times are used either alone or in combination with other phases, significant improvements in depth estimates can be obtained.In this study, we demonstrated a simple but robust procedure for estimating earthquake depths by using a stacked array beam to enhance a depth phase signal from a small scale dense seismic array at teleseismic distance is es-pecially effective for moderate-size aftershocks.
Our approach uses high quality depth phase data from a teleseismic array to avoid the tradeoff between epicenter and source depth determination.The results of this study show that most depths of the off-shore Pingtung earthquake sequence were slightly shallower than that reported by CWBSN.This discrepancy can be accounted for by the use of different earth models employed for earthquake locating.However, the later phases reflected from surface were never considered by CWBSN for earthquake locating.Though depth phases re-  corded at 85 degrees away travel the lowermost mantle (the most heterogeneous region), the downward path of pP is almost the same as P. Therefore, we do not consider a heterogeneous effect.Moreover, the layering beneath the array is horizontal and the ray paths almost perpendicular to surface so that little or no azimuthal travel time variations exist.Analysis of this study provided a new constraint for source depths evaluation of the Pingtung earthquake sequence and for further discussion of tectonic implications of this earthquake sequence.
To identify the source process of Ev2, we used numerical simulations to compare array observations (Figs.8a, b).This may support the contention that the rupture continued for a while with a major rupture occurring later than its initial onset.In spite of the initial point of the observed waveforms which indicated that the P phase was significant, the following pP phase was not clear.Therefore, we decided to choose the relatively large amplitude as a reference point to compare with synthetic waveforms and determine earthquake  depths (more detail procedure has been reported in previous section).The simulated results fit well with the observed time difference between two major P and pP phases and indicate that the source depth of Ev2 is 56 km which was deeper than that determined by the CWBSN (of about 50 km).Previous analysis of the other 5 events indicated that the source depths determined in this study are uniformly shallower than that located by CWBSN.We suggest the initial rupture depth of Ev2 was shallower than 50 km.In fact, we determined the source depth by using the pP-P time delay of relative large amplitude point, not its initial onset of earthquake fault, which should represent the depth of the rupture centroid.The discrepancy between both analyses can be reasonably interpreted as explained by the method that CWBSN uses to pick phase arrival times for its initial rupture point and obtain the initial onset depth, so the result is deservedly dissimilar with this study.Unfortunately, we do not have a horizontal location resolution using this method.If we agreed the source depth of the large asperity determined in this study, we can suppose that the rupture began at a shallow depth, continued to grow to a depth of 56 km where the largest rupture occurred and released most energy.After detailed examination of seismic data from the teleseismic array and near source strong motion seismometers (Figs. 7,8), we suggest that the rupture properties of Ev2 were compatible with our current interpretations of this earthquake sequence.

Fig. 2 .
Fig. 2. (a) Map showing the locations of the 1 st and 2 nd events, great circle path and seismic array used in this study.The focal mechanisms come from a Global CMT solution.(b) Configuration of Yellowknife Array.This array is deployed in a cross shaped configuration, the branches are oriented in W-E (R-branch) and N-S (B-branch).The short-period stations using in study are indicated by triangles.Black lines are ray paths of analyzed events approaching array stations.

Fig. 4 .
Fig. 4. A comparison between the individual station waveform record (black line) and the array beam (red line) for the 4 th event.The 1 st trace is the beam which stacked 18 stations' waveforms by the slant stack method.

Fig. 3 .
Fig.3.The left panel shows ray paths for P (in solid line), PcP (P wave reflected at the core, in dash line) and surface-reflection phases (thin solid line) from the source (black star) to receiver (black triangle).The right panel shows the detail ray paths of P and pP in source side where the black star represent source.The downward reflected ray paths, named depth phases, are nearly parallel to the P phases and the mainly difference on path is between surface and hypocenter.

Fig. 5 .
Fig. 5. (a) Comparison of the array beam (red line) for Ev4 (26 December 2006, 15:41 UTC) in Table1and the synthetic waveforms calculating different depths (black line, depths from 20 to 90 km).The time difference between the main arrival P and pP is best fit by the 8 th trace synthetic seismogram with a depth of 34 km.(b) The color image shows a coherence map that presents the consistency of the time difference of P and pP phases within observed and calculated seismograms.The maximum value (marked by a black cross) on the map means the minimum residual between the pP-P time differences of observed and synthetic records which corresponds to the optimum depth is about 34 km.

Fig. 7 .
Fig. 7.The upper map shows strong motion records of the four nearest stations (inner circle of the lower map).T1 and T2 are theoretical predictions of P and S wave arrival times, respectively.All of the records clearly present direct P arrivals.The S waves can be recognized in expectant times but with complex vibration following.The red dashed line marked the time of S with maximum amplitude.The lower map shows the locations of the strong motion stations and epicenter of event Ev2 (Large red dot).

Fig. 8 .
Fig. 8. (a) The 1 st blue trace is the stacking beam of the recorded YKA waveforms for Ev2; the 2 nd red trace is the array beam by low-pass filter with a signal less than 0.5 Hz; the other black traces are synthetic waveforms at different depths from 26 to 90 km with the same frequency band of the filtered trace.The dash line area has its zoom-in plot shown on the insertion picture.Moreover, the low-pass beam plot is also overlapped on the inserted plot to have a clear comparison.(b) The color image shows the consistency of the time difference of P and pP phases within observed and calculated seismograms.The estimated optimum depth is about 56 km.