Slip Partition of the 26 December 2006 Pingtung , Taiwan ( M 6 . 9 , M 6 . 8 ) Earthquake Doublet Determined from Teleseismic Waveforms

Teleseicmic vertical P waves of an earthquake doublet occurring on 26 December 2006 offshore southwestern Taiwan are analyzed to recover spatial slip distribution. Ground-displacement records for period 0.01 0.5 Hz are inverted using linear inversion for finite fault analysis. Waveform inversion results show that a dual segment fault of different strike is necessary for the first event to explain teleseismic waveforms. The trend in fault geometry, thus, inferred was coincidentally consistent with the bathymetry in the fault rupture area, suggesting the possible tectonic structure controlling fault propagation. Both earthquakes show larger slip in the asperities near the hypocenter and there are two main asperities of similar moments. The first event ruptured along a high angle eastern dipping plane while the second event ruptured along a low-angle western dipping fault plane. The spatial distribution of the Pingtung earthquake doublet shows a compensative and conjugate relationship, and the seismic moments of the earthquake doublet were 2.7 10 and 1.8 10 dyne-cm, respectively. The stress drops corresponding to the asperities of the two earthquakes have values of 3 to 4 bars, individually, similar to the values observed in other subduction zones.


INTRODUCTION
Two earthquakes of similar local magnitude (7.0) occurred on 26 December 2006 offshore southwestern Taiwan near the Pingtung area.These two earthquakes occurred about 8 minutes apart.One occurred at 12:26:21 UT, and the other at 12:34:15 UT.We consider these two earthquakes of similar magnitude and location as a doublet event.Hereafter, we refer to this doublet as the 1226 event for the earthquake occurring at 12:26, and the 1234 event for that which occurred at 12:34.Locations determined by the Central Weather Bureau (CWB) show these two events were very close to each other (see Fig. 1a).Relocated aftershocks determined by Wu et al. (2008) using three-dimensional velocity structure for the time period 26 December 2006 to 31 January 2007 were distributed primarily around the epicenter of the 1234 event and exhibited a NW-SE trend in distribution.
The Pingtung earthquake doublet is interesting both for its occurrence and location offshore of southwestern Taiwan, which is an area that is relatively less active compared to northeastern Taiwan.This doublet occurred in a region where the Eurasia Plate (EUP) subducts eastward beneath the Philippine Sea Plate (PSP) in a subduction zone along the Manila Trench (MT). Figure 1a shows regional background seismicity and focal mechanisms, including the doublet.Distribution of background seismicity and focal mechanisms along profiles a1 (Fig. 1b) and a2 (Fig. 1c) are also shown.The focal mechanisms shown above are the CMT solutions determined from the Broadband Array in Taiwan for Seismology (BATS).The focal mechanisms of the doublet, determined by BATS, show that the 1226 event is a normal fault with right-lateral component and the later one, the 1234 event, is a left-lateral mechanism.The focal mechanisms determined by other agencies are also shown in Fig. 1a for reference.These generally agree on a normal faulting mechanism for the 1226 event, and a strike-slip mechanism for the 1234 event.Background seismicity is shown for events with local magnitude greater than 2.5 for the years 1973 to 2006.From a historical perspective of background seismicity, seismic activity in southwestern Taiwan is less active.A subduciton pattern has been delineated from profiles a1 (Fig. 1b) and a2 (Fig. 1c).The Pingtung earthquake doublet is located at the front of the subduciton zone in the outer-rise region.These M7.0 doublets focus attention on the potential for seismic activity in the Manila Trench subduction zone.To understand the character of the source process and its possible tectonic setting, we investigate the spatial slip distribution of the source rupture and the possible association of these two earthquakes.In addition, the source parameters obtained from the results of the finite-fault source model can be of benefit to future seismological and engineering studies.
In this study, we use teleseismic data for inversion of the spatial slip distributions of the Pingtung earthquake doublet.The extent of azimuthal coverage of teleseismic stations to this doublet determines the robustness of inversion results.Although strong motion data might provide higher frequency information, the limited azimuthal coverage of strong motion stations and the complex subduction zone wedge structure of this offshore event create difficulties in accurately calculating Green's functions.However, by inverting teleseismic P waves, the finite-fault inversion scheme has been successfully used to recover slip distribution on the fault plane in previous studies (Heaton 1982;Hartzell and Heaton 1983;Mendoza and Harzell 1988;Medoza 1995;Lee and Ma 2000).Thus, the slip distributions of both earthquakes were modeled using the telesismic waveform data of P-waves to obtain the kinematic rupture pattern and interaction of the Pingtung earthquake doublet.

DATA, FOCAL MECHANISM, AND METHOD
Teleseicmic vertical P waves were used to invert the slip pattern of the Pingtung earthquake doublet.The teleseismic stations of the 1226 and 1234 events used in this study and corresponding parameters are listed in the Table 1 and  Table 2, respectively.Due to the short time interval between the occurrence of these two earthquakes, S waves of the 1234 event were coupled with the multiple phases induced by 1226 event.For comparison of these two earthquakes, we discard S waves for data consistency in this study.In addition, stations at distance from the epicenter of the 1234 event are chosen to avoid phase interference from the 1226 event.We chose records with good data quality and station coverage of the earthquakes.Additionally, we only consider stations with epicentral distances between 30 to 90 degrees to avoid structure complexities.Instrument response was deconvolved from the original recordings to acquire true ground displacement.Then, the data were bandpass filtered between 0.01 and 0.5 Hz.By considering the similarity in waveforms at these stations, a time window of 50 sec starting at 5 sec before initiation of P waves was adopted with a sampling rate of 0.25 sec for both earthquakes.Overall, 16 stations were used for the 1226 event and 15 stations for the 1234 event, respectively.The station distribution and the teleseismic waveforms of the two earthquakes examined in this study are shown in the Fig. 2. Earthquake source parameters had been resolved soon after the earthquakes by various institutions such as CWB, 570 Yen et al.BATS, the Global CMT Project, and the US Geological Survey (USGS) as listed in Table 3. Focal mechanisms from these institutions reveal similar fault types.The 1226 event is more related to normal faulting and 1234 event to strikeslip.The polarities of telesismic P waves were plotted in Figs.3a and b for lower hemisphere equal area projections of these two earthquakes, respectively.Before setting the adaptive dimension and geometry for finite fault inversion, we first try to determine the optimum focal mechanisms of the doublet.We reference the results of the various organizations mentioned above, and adjust them to obtain optimal ones through the fitting of the polarity of P waves.Finally, an optimal solution for the focal mechanism with strike (335°), dip (75°), and rake (-98°) was obtained for the 1226 event (Fig. 3a).Similarly for the 1234 event, a focal mechanism with strike (160°), dip (52°), and rake (2°) was determined as the optimal solution with good fits of P-wave polarities as shown in Fig. 3b.We used the generalized ray method to calculate pointsource responses of P waves (Langston and Helmberger 1975;Langston 1978), and the velocity structure comprised of three layers over a half space (Wu et al. 2007) was utilized for Green's function calculations (Table 4).A constant t* for the attenuation operator was adopted properly (Lay and Wallace 1995) at a value of 1.0 sec for P waves (Futterman 1962).Green's functions were calculated for a 2 sec source time function with an equilateral triangle.The Green's functions were also bandpass-filtered (0.01 to 0.5 Hz) and resampled (0.25 sec).
Based on empirical relationship of Wells and Coppersmith (1994) and the aftershock distribution, we estimate the possible fault length and width of rupture area for a planar finite fault.We divided the finite fault into several subfaults with a point source at the center of the subfault.Each subfault had dimensions of 5 ´5 km with a given fixed focal mechanism.Two nodal fault planes of the optimal solutions for the aforementioned focal mechanisms were examined.The optimal fault plane was chosen by considering the fit of synthetic outcomes to those of observations.We reference locations from the CWB as being the initial rupture points.The rupture velocity value was set as a constant.A linear least-squares fit was implemented to obtain the slips of each subfault on the fault plane.The linear inversion equation is expressed as: where A is the matrix of Green's functions, b is the ob-served data vector, and x is related to the amount of slip on each subfault.We utilize this equation in an over-determined manner to avoid instability.Nonnegative and smoothing constrains were applied following the procedure of Hartzell and Heaton (1983).For estimating reliability, a variance value of misfit is defined as (Ax-b) 2 /b 2 .A minimum misfit value is estimated to acquire the best solution for obtaining the slip on each subfault.

Spatial Slip Distribution for 1226 Event
A finite fault model of length 120 km and width 100 km was divided into 528 subfaults of dimension 5 ´5 km for the 1226 event.We consider a rupture velocity of 2.7 km sec -1 in this model.This value is about 60% of the shear-wave velocity near the hypocenter.Higher and lower values of the rupture velocities were tested; however, these values yielded results with less satisfactory waveforms and sparse spatial slip distribution.
By examining the two nodal fault planes of the optimal focal mechanism from the waveform fitting through the inversion procedure, a NW-SE strike with an eastern dipping plane was chosen.The spatial slip distribution over the optimal fault plane is shown in Fig. 4a.The slip distribution is concentrated mainly in two regions.The rupture propagated from the hypocenter northward along the strike direction.A large asperity developed near the hypocenter with a maximum slip of 1.95 m.Other sparse asperities developed in the region 30 to 90 km north of the hypocenter.Synthetic waveforms less satisfactorily resemble the later portion of observed waveforms as shown in Fig. 4b.The less satisfactory waveforms of the synthetics were found to be for stations in the  northwest azimuth.The waveforms with black dots in Fig. 4b denote the stations with poor waveform fitting between the observed and synthetic waveforms in the later 20 -40 sec.
The analysis shows that the contributions of the spatial slip distribution for the northern and deeper slips were pri-marily from the later portion of waveforms.Unsatisfactory waveform fitting may be the result of an inadequately described mechanism.In other words, a likely explanation is inadequate geometry describing the later half of the fault plane.We, then, introduce a dual segment fault model.Fault length and width of the first asperity near the hypocenter as the first segment is fixed.Different strike angles based on the focal mechanism of this first segment as an initial value were considered and examined at every 5 degrees for the 2 nd segment.The 2 nd segment constituted a length 30 km north of the hypocenter along the strike direction to 90 km.Comparing various solutions by misfit, we obtain the optimal result for the 2 nd segment as being of strike N5°E.
The slip distribution of this dual segment fault is shown in Fig. 5a.The rupture propagates unilaterally from the hypocenter to the north and then in a down-dip direction along the strike.There was no significant discrepancy in the slip pattern near the hypocenter in the first segment as shown in the initial single segment fault model.The asperity has maximum slip of 2.3 m at a depth of 35 -60 km.In the 2 nd segment, one concentrated asperity was obtained in a deeper portion of the fault at a depth of about 65 -95 km with maximum slip of about 1 m.A spatial and temporal interval exists between these two significant asperities.This feature can be observed directly from observed waveforms at stations DZM, SNZO, and TAU.They show that the waveforms were interrupted after about 15 sec for several seconds, then, other energy arrives after approximately 20 sec with broader waveforms.
Significant improvement in waveform fitting (misfit of 0.48) was found once the dual segment fault model was introduced as shown in Fig. 5b.Allowing for the later portion of the waveforms can be now explained very well.This finding reflects the importance of introducing two-segment fault geometry into this model.However, the stations in the southeastern direction (DZM, SNZO, and TAU) show relatively smaller amplitudes compared to the observations from the 2 nd asperity.This might suggest a smaller change in slip direction due to the 2 nd asperity compared to that of the first asperity.A seismic moment of 2.7 ´10 26 dyne-cm was obtained for this model.

Spatial Slip Distribution for 1234 Event
For the 1234 event, a finite fault model of length 90 km and width 80 km was divided into 288 subfaults of dimension 5 ´5 km.The hypocenter as determined by the Center Weather Bureau (CWB) was used (see Fig. 6a).A constant rupture velocity of 2.7 km sec -1 (as with the 1226 event) was considered in this calculation.The aftershock distribution (Fig. 1) suggests a west-dipping plane as being ruptured plane (Fig. 3b).We consider 15 stations for P-wave modeling to obtain the spatial slip distribution over this western dipping plane.The spatial slip distribution for the 1234 event is shown in Fig. 6a.From the slip pattern over the fault  plane, we found that the rupture propagates slightly upward and unilaterally to the northwest.A large asperity with maximum slip of about 1 m appears in the vicinity of the hypocenter and another 20 km north of the hypocenter.
A comparison of the observed and synthetic wave-forms from the dislocation model of the 1234 event is shown in Fig. 6b.Generally, synthetic waveforms explain the observation waveforms well with a misfit of 0.5.A seismic moment of 1.8 ´10 26 dyne-cm is derived from the inversion.

DISCUSSION
The occurrence of earthquake doublets has been found in different subduction regions around the world.Several studies had been carried out to understand the characteristics of doublets in terms of seismicity, stress drop, and tectonic setting (Lay and Kanamori 1980;Astiz and Kanamori 1984;Yamamoto et al. 2002;Gibowicz and Lasocki 2007).In Mexico, Middle America and the Solomon Islands, earthquake doublets have been investigated, and shown to have stress drops of about 4 to 19 bars.In the Solomon Islands, Lay and Kanamori (1980) studied three earthquake doublets including doublets with surface magnitude of more than 7.0 that occurred in 1971, 1974, and 1975.They obtained stress drops of: 17 and 19 bars for the 1971 doublet, 12 bars in both earthquakes of the 1974 doublet, and 9 bars in both earthquakes of the 1975 doublet.In Mexico, Middle America, Astiz and Kanamori (1984) showed that an earthquake doublet of moment magnitude of 6.9 in 1982 was estimated to have a value of 4 bars.
To estimate values of the stress drops for the asperities of the Pingtung earthquake doublet, we consider a circular fault model as the stress drop Ds = 7M 0 / 16a 3 where M 0 is the seismic moment; and a is the radius of a circular fault (Kanamori and Anderson 1975).Accordingly, the stress drops with respect to the two main asperities of the 1226 event were similar and estimated to be 3.19 and 3.05 bars for the shallower and deeper asperities, respectively.For the 1234 event, the stress drops with respect to the two main asperities were also similar and estimated to be 4.31 and 4.14 bars at the hypocenter-vicinity and northern asperities, respectively.The stress drops of these earthquake doublets mentioned above seem to be comparable to the values of 4 -19 bars from other subduction doublets.This seems to suggest that subduction doublet events are characterized by similar and smaller stress drops.
To understand the character of earthquake doublets in the Taiwan region, we examined the seismic catalogues of the CWB for the time period January 1973 to December 2006.We specified earthquake doublets as being a pair of earthquakes with a magnitude difference of less than 0.25 with a distance between hypocenters of less than 90 km and a time separation of less than 1 day.We delimitated events belonging to earthquake sequences following large earthquakes.Five pairs of earthquake doublets were found.The distribution and focal mechanisms of these earthquake doublets are shown in Fig. 7a; all of the doublets were located in the vicinity of two subduction zones, the southwestern Manila subduction zone and the northeastern Ryukyu subduction zone.Figures 7b and c show background seismicity and the distribution of doublets in profiles b1 and b2, which are perpendicular to the Manila and Ryukyu subduction zones, respectively.For the Manila subduction zone, in addition to the Pingtung earthquake doublet, another doublet occurred in 1973 with similar magnitude of 6.0 and a time separation of about 40 min.The focal mechanisms of this doublet were not determined.The locations of this doublet are different to the Pintung earthquake doublet as in the outer rise, but in the upper edge of the subducting slab as shown in profile b1 (Fig. 7b).The doublet event in the Ryukyu subduction zone shows a pair of earthquakes in 1990 denoted as events B and B' occurring in the region where the slab began to subduct.They both show similar thrust focal mechanisms with similar magnitude of 6.5.But, these doublets have slightly longer separation in time at about 16 hours.Events, C and C', occurred in 2001 with magnitude 6.3 and similar strike-slip foal mechanisms.Profile b2 shown in Fig. 7c indicates that this doublet might have ruptured through the overall depth regime of the subducting slab as event C occurred at a shallower depth than the deeper event C'.Another pair of events, D and D', are the well known Ilan doublet.These relatively recent events occurred in 2005 and were only 2 mins apart and exhibited similar magnitude (6.0), focal mechanisms and locations.Comparing these doublet events, the 2006 Pingtung earthquake doublet is the largest in magnitude, while the 2005 Ilan doublet had the most similar focal mechanisms and locations.We also estimate stress drops for the Ilan doublet using the spatial slip distribution determined by Tsao (2006).We obtained stress drops of 18 and 11 bars, respectively, for events D and D'.Again the values of stress drops for this doublet event were between 4 -19 bars as observed elsewhere.In view of earthquake doublets examined elsewhere (Kagan and Jackson 1999), earthquake doublets seem to occur frequently in subduction zones and exhibit similar stress drops.
We further discuss the interaction of the spatial slip distribution of the Pingtung earthquake doublet.As shown above, fault geometry of 1226 event cannot be simply explained by a single fault segment.A 2 nd segment bending eastward is necessary to explain observed waveforms.Coincidently, the bending strike of the 2 nd fault segment correlates well with the bathymetry of the earthquake rupture area (Fig. 8).This seems to suggest the controlling tectonic structure to the earthquake fault rupture.Additionally, we found that the spatial slip distribution of the 1226 event was in a region perpendicular to the Manila Trench where the Eurasian Plate subducts the Philippine Sea Plate; whereas, the spatial slip distribution of the 1234 event was in a region parallel to the Manila Trench.Similar features were also observed for the Solomon Island doublet and Mexico doublet (Lay and Kanamori 1980;Astiz and Kanamori 1984).This indicates that stress heterogeneities might exist at the regional tectonic structure between the locations of the 1226 and 1234 events.
We compare the locations of the slip asperities projected to the surface for the 1226 and the 1234 events, respectively.The relative pattern is depicted schematically in Fig. 8.The sequence of the occurrence of asperities in time is denoted from A to D. The slip patterns projected over the surface reveal no superimposing of the asperities over each other.This suggests that the stress release of each asperity was independent in time and space.But, asperities of the 1234 event were compensated to the spatial slip distribution of the 1226 event, suggesting possible energy interaction of the two events.The slip asperities of the earthquake doublet projected to the surface were also compared to the onemonth aftershock distribution, which was relocated using three-dimensional velocity structure (Wu et al. 2008).Except for asperity C of the 1234 event, which was filled with relatively numerous aftershocks, most aftershocks occurred in areas surrounding asperities.In terms of stress release, most stress was released in the main asperities, leaving little or no stress to be released during the aftershocks.Lin et al. (2008) calculated the static stress changes of the first event, and show that the occurrence of the 2 nd event was located in a region where significant static stress increased, suggesting possible stress interaction of the earthquake doublet.
Strong motion data often provides the constraints necessary for generating higher resolution in slip pattern.However, our attempt to consider strong motion data within a 576 Yen et al.(a) (b) (c) frequency band higher than 0.1 Hz failed due to limited station coverage for the offshore events and complicated 3D velocity structure in the subduction wedge.Thus, our spatial slip distributions for the doublet were mainly for the results of a seismic period larger than 10 sec waves.Lee et al. (2008) using strong-motion waveforms for spatial slip distribution were also limited in their examination of the Pingtung doublet to a frequency band of 10 sec to 100 sec similar to the frequency band used in our teleseismic waveform inversion.However, their results did not present the NE bending of the 2 nd segment in the case of 1226 event.
This discrepancy might result from not including the NW stations in the strong motion dataset of the 1226 event, as these stations to the NW azimuth in teleseimic distance were the ones that gave the constraints resulting in the requirement of a 2 nd NE bending fault.Recent studies of earthquake scaling have drawn attention to the usage of different datasets and frequency bands in affecting earthquake scaling in terms of maximum slip and fault dimension.Our studies show that a full azimuthal coverage to the event is necessary to depict more complete fault geometry, and the amount of slip is related to the frequency band of the data used.Inadequacies might become apparent in the generation of earthquake scaling relationships for different data types and different frequency bands.

CONCLUSIONS
We investigated the Pingtung earthquake doublet via teleseismic P waveform inversion to obtain the spatial slip distribution and corresponding stress drops of slip asperities.Our results show that for the 1226 event, a two-segment fault geometry is necessary to explain the observed waveform, especially for stations in the northwestern azimuthal direction.This 2 nd bending fault segment would not be evident if analysis was done only through inland strong motion stations that have no northwestern coverage of these events.The results suggest the importance of good azimuthal coverage by stations to an event for accurate fault geometry determination.The spatial slip distributions of this doublet show maximum slip at the asperities near the hypocenter with maximum slip being 2.3 and 1.2 m, respectively.For the 1226 event, the dimensions of the asperities are 45 ´40 km in the southern portion of the fault segment and 40 ´50 km in the northeastern fault segment.The seismic moment of this earthquake is 2.7 ´10 26 dyne-cm.For the 1234 event, the dimensions of the two apparent asperities are 35 ´35 km in the southeast portion and 30 ´40 km in the northwest portion of the fault.The inferred seismic moment of this earthquake is 1.8 ´10 26 dyne-cm.These results suggest that both events in the Pingtung earthquake doublet have two asperities of similar seismic moment for respective magnitudes of 6.9 and 6.8.If we combine the seismic moment of the doublet, the total energy of those asperities is similar to that of an M w 7.0 earthquake with total rupture length of 80 to 90 km.
The asperities have stress drops of 3 ~4 bars, similar to stress drops from doublets in other subduction zones.The propagating direction of the 1226 event was northward and down-dip direction was along strike.The 1234 event had a similar propagating direction toward the north along a fault plane with a southeast-northwest strike direction, almost parallel to the first segment strike of the 1226 event.The asperities of the two events do not overlap in space and are located on two distinct fault planes with one dipping eastward for the 1226 event and the other westward for 1234 event.
There is a suggestion that the spatial distribution of the Pingtung earthquake doublet shows a compensative and conjugate relationship.

Fig. 1 .
Fig. 1.(a) The locations of events 1226 and 1234 as determined by the CWB. 1 and 2 denote 1226 and 1234, respectively, shown as asterisks.Black dots indicate aftershocks from December 26 to January 31.Focal mechanisms from different agencies for the doublets are shown, including those from BATS.Inserted figure is our research region, where EUP subducts eastward beneath the PSP along the Manila Trench.Corresponding seismicity of 25 km width of (b) profile a1 and (c) profile a2.Gray dots are background seismicity for 1973 to 2006.
waveforms were acquired from the Incorporated Research Institutions for Seismology (IRIS).

Fig. 2 .
Fig. 2. Station distribution and the teleseismic vertical P waveform used for (a) 1226 event and (b) 1234 event.The asterisks and triangles denote the epicenter and location of stations, respectively.Maximum amplitude in micrometers and the azimuth in degrees of each station are also shown.

Fig
Fig. 3. P-wave first motions of (a) 1226 event and (b) 1234 event.Solid and open circles show compression and dilatation, respectively.Dashed lines display optimal fault planes for spatial slip distribution determination.Projection is a lower hemisphere equal area projection.

Fig. 4 .
Fig. 4. (a) Spatial slip distribution of the 1226 event for a single segment fault (inserted figure).The hypocenter is shown by asterisk.Slip is contoured with interval of 0.4 m.(b) Comparison of observed (black solid line) and synthetic (red dot line) waveforms for the 1226 event.Numbers at the beginning of each record indicate peak displacement amplitude of observed waveforms in micrometers.The station codes and azimuth are shown at the end of each record.Black dots indicate portion of waveforms with less satisfaction.

Fig. 5 .
Fig. 5. (a) Spatial slip distribution of the 1226 event for a dual segment fault (inserted figure).Hypocenter is shown by asterisk.Slip is contoured with interval of 0.4 m.Dotted line indicates the dividing line of the two segments.Strike direction of each segment is shown above.(b) Comparison of observed (black solid line) and synthetic (red dot line) waveforms for the dual segment fault model of 1226 event.Numbers at the beginning of each record indicate peak displacement amplitude of the observed waveforms in micrometers.Station codes and azimuth are shown at the end of each record.

Fig. 6 .
Fig. 6.(a) Spatial slip distribution of the 1234 event.Hypocenter is shown by asterisk.Slip is contoured with interval of 0.4 m.(b) Comparison of observed and synthetic waveforms for 1234 event.Numbers at the beginning of each record indicate peak displacement amplitude of observed waveforms in micrometers.Station codes and azimuth are shown at the end of each record.

Fig. 7 .
Fig. 7. (a) The distribution of earthquake pairs defined as doublets for years 1973 to 2006 in Taiwan region.Date in yyyymmddhhmm format and local magnitude in parentheses are indicated for those pairs.Five doublets as A and A', B and B', C and C', D and D', and E and E' are shown by asterisks.The profiles (b) b1 and (c) b2 are depicted with background seismicity, as shown in gray dots, of magnitude greater than 3.0.