Source Mechanisms and Rupture Processes of the 26 December 2006 Pingtung Earthquake Doublet as Determined from the Regional Seismic Records

The 26 December 2006 Pingtung earthquake doublet (both about Mw 7.0) were the two largest events with magnitudes close to 7.0, following the 1999 Chi-Chi earthquake. The locations and focal mechanisms of the Pingtung earthquake sequence provide good opportunities for scientists to investigate plate movement and deep tectonic structure in southern offshore of Taiwan. However, precise hypocenters and source mechanisms of this sequence are not definite due to poor station azimuthal coverage and a lack of information from historical seismicity in this area. Here we use regional ground motion full waveform records to investigate source rupture processes and determine fault planes of the Pingtung earthquake doublet. The finite fault inversion result of the first event shows two large asperities occurred near the epicenter. The rupture mainly originated from the south then propagated northward on an eastern dipping fault plane. The second earthquake is a bilateral rupture along NW-SE direction having a large asperity close to the epicenter which released most of its seismic energy within the first 15 seconds. The slips of the two earthquakes show compensative patterns when mapped on the earth’s surface. We show that the Pingtung earthquake doublet are large normal and strike-slip faulting which dips toward the east and west, respectively. The mechanisms determined from inversion results infer that the first event might relate to the bending of the South China Sea Plate during its subduction to the Philippine Sea Plate. The second event could happen at a pre-existing weak zone that was triggered by transferred stress after the first event.


INTRODUCTION
On 26 December 2006 (12:26:21.0UTC), a M L 6.8 earthquake occurred offshore southern Taiwan.About eight minutes later (12:34:15.0UTC), another earthquake (M L 6.9) occurred in a similar area.The hypocenters reported by the Central Weather Bureau of Taiwan (CWB) were 120.56°E, 21.89°N (depth 21.9 km) and 120.39°E, 21.95°N (depth 47.03 km) respectively.A few days later, the magnitudes of the two earthquakes were corrected to M L 7.0.The epicenter of the first earthquake was relocated to 120.56°E, 21.69°N (depth 44.1 km), south of the original epicenter, with a much deeper focal depth.These two earthquakes were the largest doublet events occurred offshore southwest Taiwan over the last century (Fig. 1).
The fault planes derived from BATS (Broadband Array in Taiwan for Seismology) CMT solution indicate normal fault striking 349°, dipping 53°, and raking -54°(with strike 119°, dip 50°, rake -128°on an auxiliary plane) for the first earthquake; and strike-slip fault striking 144°, dipping 26°, and raking -12°(with strike 119°, dip 50°, rake -128°on an auxiliary plane) for the second event.It is a challenge to look for the true rupture planes because the precise hypocenters and focal mechanisms of these two offshore earthquakes are difficult to be determined from inland records.Although the earthquake reports from US Geological Survey (USGS) and Harvard University also provided CMT solutions of the doublet events, the source information are still ambiguous due to limited resolution of the teleseismic data used.Furthermore, compared to the 1999 Chi-Chi earthquake which had nota-ble surface break (Chen et al. 2001) and well known subsurface structures (Wang et al. 2000) that helped to constrain the fault plane (Lee and Ma 2000), seismologists and geologists encounter more challenges to determine rupture planes of the Pingtung doublet events.
In this study, we try to resolve fault planes as well as rupture processes of the earthquake doublet by considering finite fault source inversion based upon waveforms from ground motion data.We apply a parallel inversion technique (Lee et al. 2006) to invert for slip distributions and rupture time histories on the possible fault planes.The inversion results can help clarify source mechanisms and rupture processes of the 2006 Pingtung offshore earthquake doublet.

DATA AND METHOD
Regional ground motion waveforms are considered in the inversion study, including data recorded by TSMIP (Taiwan Strong Motion Instrumentation Program, Wu et al. 2008) and BATS stations during the Pingtung earthquakes.Figure 2 shows the distribution of used stations and their recorded waveforms.Because parts of the near field broadband records are saturated, such as station TWKB and WLCB, we consider TSMIP records to fill in the lack of near field ground motion information.As can be seen in Fig. 2,  the stations close to the source area are mostly supplied by TSMIP.The strong motion acceleration waveforms are integrated to velocity first and then band-pass filtered between periods of 10 to 100 seconds (0.01 -0.1 Hz), followed by base-line correction and resampled at 0.5 seconds.The broadband velocity data employs the same data processing but skips the integration step.Even though we have considered most of the regional records, but as shown in Fig. 2, the stations are mostly distributed on one side of the Pingtung earthquake, i.e., on the island of Taiwan where the azimuthal coverage of stations is poor.However, by considering full waveform in the source inversion study, not only the travel times but also the radiation pattern and amplitude of all phases, can be used to constrain the inversion result.
We set the possible fault planes at depths between 1 and 80 km and the faulting areas are bounded between 20.9°N to 22.6°N and 119.5°E to121.2°E.The fault planes are then divided into numerous subfaults (5 ´5 km) by a finite fault approach.The ground motion record at a given station can be represented by a linear sum of slips weighting from these subfaults.The synthetic Green's function on each subfault to a given station was generated by the F-K method (Wang and Herrmann 1980).The velocity model used to compute the Green's functions is a 1-D Taiwan average crust model (Chen and Shin 1998).The observed data and Green's functions can be formed as a system of linear equations, where A is the matrix of synthetics, b is the data vector, and x is the solution matrix of the subfault slips.This equation is solved by linear least squares, but sometimes the solution is unstable, so the problem may be stabilized by adding linear stability constraints, including smoothing, damping, and moment minimization.The waveform misfit is defined as (Ax -b) 2 / b 2 .The smaller value of misfit indicates the less discrepancy between synthetic and observed waveforms.When the misfit is zero that means the synthetic and observation are exactly the same.This inversion method was successfully used in analyzing the slip distribution of many large earthquakes happened in the world, such as 1983 Imperial Valley earthquake (Hartzell and Heaton 1983), 1992 Landers earthquake (Wald and Heaton 1994), and 1994 Northridge earthquake (Wald et al. 1996), and 1999 Chi-Chi earthquake (Lee and Ma 2000;Lee et al. 2006;Ma et al. 2001).For large earthquakes, source time functions are usually complex within the fault.We consider the multiple-time window which allowed each subfault to slip in any one of 2.0 s time windows following the passage of the rupture front.In order to resolve the detailed variation in the waveform, each window is set to have an overlap of 1.0 s.We also assume a bit faster initial rupture velocity of 3.0 km s -1 .After the inversion by considering a multiple-time window, the rupture velocity could be varied between 0 -3.0 km s -1 .
However, an increase in the number of time windows would result in a large expansion of matrix A. Lee et al. (2006) developed a parallel nonnegative least squares (NNLS) inversion technique which essentially decomposes matrix A into different computing nodes, and thus improving the number of time windows and promoting program performance.Figure 3 is a diagram showing the multiple-time window approach and its relationship to the parallel NNLS inversion.The optimum capacity of our computer facility, using a parallel NNLS program, can solve at least 40 time windows under 10 nodes 20 CPUs IBM Cluster (IBM eServer in the Institute of Earth Sciences, Academia Sinica, IESAS).
We follow two steps to determine the weighting values of constraints used in the inversion.First, a trade-off test between misfit and seismic moment are performed to decide the weighting of minimized solution constraint.A spatial-smoothing weighting test is then applied to confirm the final slip distribution as well as to compromise reciprocal requirements between model resolution and estimation errors.Figure 4 shows an example of this procedure.In Fig. 4a, the larger minimized solution weighting reduces the amount of seismic moment but increases the waveform misfit.We obtain an optimum weighting for minimized solution constraint of 0.12.For spatial-smoothing weighting value of 0.2, the final solution is then decided (Fig. 4b).It is notable that the spatial-smoothing constraint has less influence on waveform misfit and final seismic moment compare to minimized solution constraint in this study.

The First Event (2006/12/26 12:26)
We analyzed inversion results to estimate the true fault plane.The inverted slip distributions on possible fault planes are shown in Fig. 5.The left panels are two faults determined from BATS CMT (Fault 1, Fault 2), and the right panels are fault planes from Harvard CMT (Fault 3, Fault 4).Faults dipping toward the east are shown in the upper diagrams and the lower diagrams show fault planes dipping toward the west.The relevant stations, records and inversion parameters, including the initial rupture velocity (3.0 km s -1 ), weighting of slip smoothing, damping and moment constraint, are identical in this fault plane examination.
First we analyzed the fault dips toward the east as determined from BATS CMT solution (Fault 1 in Fig. 5a).The slips on the fault plane are concentrated in several areas and the slip vectors are smoothly varied.For the auxiliary fault plane which dips westward (Fault 2 in Fig. 5b), the distribution of slip is also concentrated and part of the slip occurred beneath southern Taiwan.The slips on Fault 3 (Fig. 5c) show a similar distribution compared to the result of Fault 1 because these two fault planes are alike and dip toward the east with similar strike and dipping angle.However, for the fault plane dips toward the west determined from Harvard CMT (Fault 4, Fig. 5d), the slip shows a complex split pattern which occurs along two different depths at shallow part of the fault.In general, faults dipping toward the east have smaller misfits within four candidates.Furthermore, the slip patterns are more comparable with the distribution of aftershocks for the east dipping fault planes.It is notable that the major slip area on Fault 3 is not identical to centroid as might be expected from the Harvard CMT solution.From the inversion results and previous analyses of all four candidates, the preferred fault plane of the first event is Fault 1 which is determined from BATS CMT solution that dips toward the east.
We further investigated the best-fit slip distribution on Fault 1 by considering different weighting values in the in-version.After several tests, the final inversion result is shown in Fig. 6.Most of the rupture areas occur along north-south direction and slip toward a down-dip direction (normal faulting with a rake angle close to 90 degrees).Two big asperities are found on the fault plane, one located close to the hypocenter with a larger slip area, and another to the south of the hypocenter which has the largest slip amplitudes of about 319 cm.The slips on northern asperity are mostly deeper than 35 km; however, the slips on the southern asperity are shallower and are about 25 to 35 km.It is notable that rupture in the northern part of the fault becomes smaller and their slip vectors have more strike-slip component which slips toward the north.The observation and synthetic waveforms have good fits in general except the stations close to the source area and located in the Western Plain (such as sta-   tion KAU038, KAU042 and TAIB).The total misfit of waveform is 0.589.Total moment determined from the inversion result is 4.26 ´10 26 dyne-cm which is equal to about M w 7.02.

The Second Event (2006/12/26 12:34)
The focal mechanisms of the second event determined from BATS and Harvard CMT are similar with one fault plane dipping westward and another one approaching a vertical strike-slip.We examined the slip distribution of these four possible fault planes as shown in Fig. 7. Again, the left diagrams are two fault planes determined from BATS CMT (Fault 1, Fault 2), and the right diagrams are fault planes from Harvard CMT (Fault 3, Fault 4).Faults dipping westward are shown in the upper part and the lower part shows vertical strike-slip fault planes.All relevant stations, records, and inversion parameters are the same as that used in the first event.We projected slips of two vertical strike-slip faults on fault planes dipping 55 degrees in Fig. 7 for easy comparison.
The slip distributions on two vertical strike-slip fault planes are scattered (Figs.7b, d) and their slip vectors are disordered.In general, ruptures on fault plane are usually not present with this behavior.On the opposite, fault planes dip toward the west show concentrated slip patterns, especially the result determined from Harvard CMT solution (Fault 3).The slip on this fault plane is not only centralized but the slip vectors are smoothly varied.The slip pattern and aftershocks are comparable in which both follow a NW-SE direction.Furthermore, Fault 3 also comprises the smallest misfit within four candidates.From these analyses, the best fault plane of the second earthquake determined from the waveform inversion is the fault derived from Harvard CMT solution which dips toward the west.
Fault 3 is further used to investigate the final slip distribution of the second event.The result is shown in Fig. 8. Three larger asperities are found, the largest one just adjacent to the epicenter which also has the largest slip on the overall fault plane (260 cm); the other two are smaller and occur beside the first asperity.All slips are distributed along a NW-SE fault plane and slip toward strike direction.In general, the slips on all three asperities occurred below depths of 30 km.The waveform misfit between synthetic and observation is 0.564 which is better than that of the first event.However, the waveform comparisons for stations close to the source and located in the Western Plain are still not much improved.The total moment obtained from the inversion result is 3.97 ´10 26 dyne-cm (M w 7.0) which is similar to the first event.

DISCUSSION 4.1 Rupture Snapshots
The rupture time history can be reconstructed when we apply a multiple-time window technique in a source inversion study.The cumulate snapshots of the two events are shown in Fig. 9. Starting from the first event (Fig. 9a), the initial slips mostly occurred near the epicenter within the first 5 seconds.Between 5 to 10 seconds, a big asperity occurred to the south of the epicenter (at about 21.5°N) which continued to slip for about 5 seconds.The rupture behavior at this moment is in agreement with the result of initial rupture processes analysis using near source strong motion records (Huang et al. 2008).This activity resulted in the largest slip of the first event.From 15 to 25 seconds, the rupture began to propagate northward and during that time the second asperity in the middle part of the fault plane arose.After that, the rupture continued to propagate north with smaller slips between 25 to 35 seconds.The rake angles of slips varied gradually from -90 to about -30 degrees during this period.The moment rate function of the first event is shown at the right of the diagram in Fig. 9a.The first asperity occurred at about 5 to 10 seconds and the moment rate function dropped down.At about the 15 th second, the second asperity occurred and was followed by another peak moment release when the rupture propagated northward.The total rupture duration time of the first event was approximately 35 seconds that released a seismic moment of about 4.26 10 26 dyne-cm.Eight minutes later, the initial rupture of the second event occurred which was close to northern end of the first event (Fig. 9b).The slip in the first 5 seconds was quiet, but large rupture areas were suddenly present at the 10 th second.From 5 to 10 seconds, a large amount of slips arose around the epicenter.This phenomenon indicated that the second earthquake might begin with a quickly ruptured behavior.This area continued to slip until the 15 th second, resulting in two asperities, one close to the epicenter and the other northwest of the epicenter.The rupture becomes quiet until the 20 th second, and then southeast fault plane began to slip until the 35 th second which grew into a third asperity.We analyzed the moment rate function of the second event at the right of the diagram in Fig. 9b.The moment rate function is simpler to compare to the first event.The first two asperities released most of its seismic energy between 5 to 15 seconds then the moment suddenly drops down.The third asperity released energy smoothly between 20 to 35 seconds then gradually grew quiet.The total seismic moment of the second event was about 3.97 ´10 26 dyne-cm which mostly released energy within the first 15 seconds.In general, the second earthquake ruptured bilaterally while the first earthquake ruptured toward the north.

Total Slip Pattern
We projected the total slips of the two events to the ground surface, a continuous slip pattern are reserved as shown in Fig. 10.The two asperities of the first event are located at the southern part.There are almost no aftershocks located near the asperities.Although there are some slips occurred when the rupture propagated to the north, in general the rupture stopped at about 22°N.It is notable that the rupture of the second event originated near this area where the largest asperity as well as a large amount of aftershocks occurred here.Mainly, the major rupture areas of the Pingtung earthquake doublet are not overlapped with each other.On the other hand, they present a compensative relationship between asperities.This result is similar to another source inversion result determined from teleseismic waveforms (Yen et al. 2008).From south to north, the depths of these asperities become deeper and deeper which vary between 25 and 45 km.This phenomenon might relate to depth variation of the Moho southwest offshore of Taiwan.Furthermore, there are no definite reports (Dennis et al. 1995;Susan and Dennis 1995a;Susan et al. 1995b;Dennis and Susan 2005) indicating discontinuous subsurface structures close to the source area which could potentially stop a rupture, and then result in the second event.Stress trigging would be another candidate to explain the relationship between the earthquake doublet (Lin et al. 2008).The second event could happen at a pre-existing weak zone that was triggered by transferred stress after the first event.To realize these possibilities, more insight investigation should be taken in both seismological and geological studies.

Source Rupture Model Parameters
The source rupture model parameters de-termined from inversion results are provided in Table 1.Major rupture areas of the Pingtung earthquake doublet were 50 ´35 km 2 and 65 ´30 km 2 , respectively.The seismic moment of the first earthquake was slightly larger than the second event which was about M w 7.02 and 7.0 respectively.The average stress drops of the two events are similar and close to 70 bars.This result implies that the Pingtung earthquake doublet might be intraplate earthquakes as suggested by Kanamori and Anderson (1975).Although the resolution of our inversion might be limited due to a station azimuthal coverage problem, it is still notable that waveform inversion results prefer a normal faulting dips toward the east of the first event.Most of the slips are occurred at depths deeper than 35 km.A large normal fault does not occur in the subduction zone frequently.It is possible that the plate bent and broke during its subducting process (Davies and Blanckenburg 1995).This implies that the Pingtung Earthquake could be a bending-related event which broke the deeper curst and even cut into the uppermost mantle when the South China Sea plate subducts toward east to the Philippine Sea Plate.

CONCLUSIONS
We used seismic records to investigate rupture processes of the Pingtung earthquake doublet.The finite fault source inversion result can help to clarify the fault plane as well as the related tectonic process.The first event shows two large asperities, one occurred close to the epicenter and another located to the south of the epicenter.The rupture mainly originated from south to north on an eastern dipping fault plane.The second earthquake ruptured bilaterally along a NW-SE fault plane which dipped toward the west.A large asperity occurred close to the epicenter that released most of the seismic energy within the first 15 seconds.Slips of the two earthquakes show a compensative relationship when they are mapped to the ground surface.We conclude that the Pingtung earthquake doublet are large normal and strike-slip faults, dipping toward the east and west, respectively.Inversion results indicate that the first event might relate to the bending of the South China Sea Plate during its subduction to the Philippine Sea Plate.The second event could happen at a pre-existing weak zone that was triggered by transferred stress following the first event.

Fig. 1 .
Fig. 1.The locations of the Pingtung earthquake doublet.Hypocenters and focal mechanisms reported by different institutions are shown.Stars indicate the epicenters of the earthquake doublet reported by CWB.The beach-balls labeled by 1226 indicate the first event, and 1234 indicate the second event.Depth information of hypocenters are shown in brackets.Aftershocks reported by CWB which occurred within 5 days after the mainshock are plotted according to the magnitude as indicated in the legend at the right.

Fig. 2 .
Fig. 2. Stations and records used in this study.(a) First event (12:26); and (b) Second event (12:34).The TSMIP and BATS stations are indicated by squares and triangles, respectively.The East, North and Vertical components, velocity waveforms are shown.All the records are band-pass filtered between 10 and 100 seconds.Stars show the locations of the Pingtung earthquake doublet.

Fig. 3 .
Fig. 3.A Multiple-time window approach and its relationship to the parallel NNLS inversion.(a) A schematic diagram illustrating the use of multiple-time window.(b) A multiple-time window with matrix decomposed in the Parallel NNLS inversion.Message Passing Interface (MPI) is applied as the communicant between the computing nodes in the parallel computing process.
Fig. 4. Trade-off curves of the misfit and seismic moment for various weighting values: (a) minimized solution constraint, (b) spatial-smoothing constraint.

Fig. 5 .
Fig. 5. Inversion tests on four possible fault planes of the first event.(a) Fault 1, the BATS CMT solution dips toward the east; (b) Fault 2, the BATS CMT solution dips to the west; (c) Fault 3, the Harvard CMT solution dips toward the east, and (d) Fault 4, the Harvard CMT solution dips westward.The focal mechanisms are indicated by beach balls; strike, dip and rake angles are labeled in brackets.Misfits between observation and synthetic waveforms are shown at the lower right accordingly.The depths on the fault planes are labeled by blue texts.The slip values, slip vectors and aftershocks are described in the lower legends.

Fig. 6 .
Fig. 6.Final inversion results of the first event.The slip values are shown using a color pattern as indicated at the right.Data (black) and three-component synthetic waveforms (color lines) are plotted on the map of the south Taiwan area.The focal mechanism of the main shock is indicated by the beach ball.Both the data and the synthetic seismograms are band-pass filtered between 10 and 100 seconds.

Fig. 7 .
Fig. 7. Inversion tests for four possible fault planes of the second event.(a) Fault 1, the BATS CMT solution dips toward the west; (b) Fault 2, the BATS CMT solution shows a vertical dip strike-slip faulting; (c) Fault 3, the Harvard CMT solution dips toward the west, and (d) Fault 4, the Harvard CMT solution showing a vertical dip, strike-slip faulting.The slips of the two vertical strike-slip faults are projected to fault planes dipping to 55 degrees in Figs.7b and d for easy comparison.

Fig. 8 .
Fig. 8. Final inversion results of the second event.The descriptions are the same as that in Fig. 6.

Fig. 10 .
Fig. 10.Projected slip distribution of the Pingtung earthquake doublet.The epicenters and focal mechanisms are indicated by stars and beach balls, respectively.The slips of the two events are shown by contours with different colors as described in the legend at the right.Aftershocks are also plotted in the figure for comparison.