Precursory Seismic Activation of the Pingtung ( Taiwan ) Offshore Doublet Earthquakes on 26 December 2006 : A Pattern Informatics Analysis

We report the retrospective analyses of frequency-magnitude distributions (FMD) and pattern informatics (PI) for the Pingtung offshore doublet earthquakes having magnitudes 6.7 and 6.4. The frequency-magnitude distributions of earthquakes that occurred in the study area preceding the Pingtung doublets demonstrate the activation of moderate events since March 2004, about three years before the Pingtung doublets. In addition, in the calculated PI map, the hypocentral area of the doublets exhibits the signatures of anomalous activity associated with the precursory seismic activation. Our study therefore suggests that, by inspecting temporal variations in the frequency-magnitude distributions of earthquakes, together with the hotspot locations with intense seismicity changes in the PI map, it is possible to shed light on the relevant preparation processes of the forthcoming large earthquakes.


INTRODUCTION
Taiwan is located in the interaction zone of the Eurasia and Philippine Sea plates.The interaction of two plates generates many earthquakes.Some of these major earthquakes, such as the Pingtung offshore doublet earthquakes (Fig. 1a) that struck southwest Taiwan on 26 December 2006, caused significant damage and injuries.The Pingtung offshore doublet earthquakes (hereafter called the Pingtung doublets) not only caused damage and injuries on land, but also damaged several undersea fiber-optic cables used to route internet traffic and telephone services between Taiwan, Hong Kong, Japan, China, South Korea, Philippines, Malaysia, Singapore, and Thailand (Weng et al. 2007).The first main shock occurred at 8:26 p.m. on 26 December 2006 (local time) with a normal faulting mechanism having local magnitude 6.7.About 8 minutes later, the second main shock occurred at 8:34 p.m., with a strike-slip mechanism having magnitude 6.4.The Pingtung doublets were located in a historically inactive region where no earthquake with magnitude larger than five had occurred since 1973.
Understanding the mechanism of the occurrence for the Pingtung doublets may help us to better realize the tectonics for the region offshore from southwest Taiwan.But in addition, mitigating seismic risk demands that seismologists attempt to forecast or predict when and where large earthquakes will happen.Unfortunately, the initiation process of earthquakes is complex and remains basically unsolved, raising skepticism among the seismological community for the possibility of predicting future earthquakes from observed seismicity.Nevertheless, by the use of sophisticated statistical analysis, precursory phenomena related to seismic quiescence and activation for many major earthquakes have been frequently identified and reported (Wyss and Habermann 1988;Bowman et al. 1998;Jaume and Sykes 1999;Zoller et al. 2002;Chen et al. 2005;Nanjo et al. 2006).
One of these statistical analysis methods, Pattern Informatics (PI) (Rundle et al. 2003) was developed to forecast the possible locations of forthcoming large earthquakes by investigating the space-time patterns of past earthquakes.Precursory patterns obtained from the PI method are typi-cally associated with both seismic quiescence and activation.As revealed by the PI method, dramatic changes in seismicity rate during the precursory time will be appear as a characteristic anomalous pattern related to the forthcoming events.The precursory time can be identified by analyzing the frequency-magnitude distributions (FMD) of earthquakes (Rundle et al. 2000c;Chen 2003).In this study we combine the FMD analysis with the PI method to reveal the precursory phenomena associated with the Pingtung doublets.

DATA
The earthquake catalogue published by the Central Weather Bureau (CWB) of Taiwan was used in this study.The CWB began to install a new local seismic network, CWBSN, at the beginning of the 1990s.For high performance seismic monitoring, digitally telemetered and recorded short-period velocity sensors as well as force balanced accelerometers were installed at all 73 stations around the Taiwan area.The CWB finished the set-up jobs for the network at the end of 1993.Since 1994, the monitoring operation has been changed to a continuously recording mode along with the manual identification of events, thus greatly enhancing the detection sensitivity for earthquakes.We therefore ignore the earlier catalogue data before 1994.Furthermore, by fitting the FMD of the earthquake catalogue (Fig. 1b), the catalog completeness magnitude M c for the study area can be conservatively determined as M c = 3.2.Also shown in Fig. 1c is the depth distribution of earthquakes occurring in the study area, with magnitude larger than M c .About 80% of these earthquakes occurred at depths shallower than 30 km and there exists an abrupt decrease in the frequency of earthquake occurrence at a depth of 30 km.Therefore we have divided the catalogue into two groups, those shallower and those deeper than 30 km.The Pingtung doublets, which are in the latter group, occurred at the depths of 44.1 and 50.2 km respectively.The depth range from an upper bound of 30 km to a lower bound of 100 km was consequently considered in the following PI analysis.We discuss below the effect of depth selection on the PI map.

MODIFIED PATTERN INFORMATICS METHOD
The PI method was inspired by the concept of pattern dynamics (Mori and Kuramoto 1997;Rundle et al. 2000a, b).Though the real dynamics D t of a complex system is unknown, the apparent dynamics PD t of space time patterns of events can be realized by observable variables (Fig. 2) (Tiampo et al. 2002b).In a complex system, both the real underlying dynamics D t and the real system state variables s (x, t) control the physics of the system and produce the observable pattern state variables y (x, t).Since the pattern variables y (x, t) and y (x, t + Dt) are both observable it  becomes possible to construct an approximate dynamics, the pattern dynamics, that associates a future pattern state y (x, t + Dt) with an earlier state y (x, t).The operator PD t of the pattern dynamics constructed in this way is assumed to be linear over the relatively small, observed time interval Dt (Rundle et al. 2000a).
Pattern informatics is a promising method to detect and display the characteristic precursory patterns before threshold events such as earthquakes.Since strong space-time correlations represented by seismicity are responsible for the cooperative behavior in the earthquake fault system, earthquake seismicity can therefore be used to represent the state vector in a phase dynamical system (Rundle et al. 2000a;Tiampo et al. 2002a).Tiampo et al. (2002b) suggested that seismicity can be described by pure phase dynamics (Mori and Kuramoto 1997;Rundle et al. 2000a, b), in which the important changes in seismicity are mainly associated with rotations of the state vector in a high-dimensional Hilbert space (Fukunaga 1970;Holmes et al. 1996;Mori and Kuramoto 1997;Rundle et al. 2000a).Accordingly, temporal variation in seismicity corresponds to a time-dependent drift in the phase angle over a pre-defined period, i.e., the "change interval."In the PI calculation, the vector difference over the change interval can be derived from the mean drift angle, and then the probability amplitudes of threshold events in the future can be represented by the PI value.
The procedures for the PI method (Chen et al. 2006a) are described in the following: 1.The study area is divided into many grids with the dimension of 0.1°´0.1°. 2. Four time parameters are defined.First, t 0 denotes the beginning time of the catalogue.As we mentioned above, t 0 is chosen as 1994/1/1 in this study.Second, t 1 and t 2 mark the start and end of the change interval, respectively.In general t 2 is identified as the current time.
In our retrospective PI calculation of the Pingtung doublets we have chosen 2006/5/1 (about six months before the Pingtung doublets) as t 2 .As to the start of the change interval, we explain the determination of t 1 with the FMD analysis in the following section.Finally, a sampling reference time t b shifts between t 0 and t 1 , which is used for eliminating the strong effect due to a large number of clustering aftershocks occurring at some sites during the change interval and the random fluctuations due to noise that occur as well.3. The seismic intensity is defined as the average number of earthquakes with magnitude larger than M c that occurred in a grid box x i and its eight neighboring boxes (the Moore neighbors): (1) where n(x i , t) is the number of events occurring at location x i (or its Moore neighbors) and time t.Thus the intensity change during the change interval, i.e., DI , indicating the relative possibility of large threshold events, is computed and normalized to the maximum value of P (x i ). 6.In a PI map the mean-subtracted index of occurrence probability, i.e., DP(x i ) = P(x i ) -m P , is then color coded and plotted.Here m P is the mean of P(x i ) over all boxes and can be considered to be the average background probability for the whole area.

RESULTS
According to the self-organized spinodal (SOS) hypoth- For the detailed explanation please refer to Section 3.
esis proposed by Rundle (2000c) the earthquake cycle from one characteristic earthquake to another can be divided into three periods as shown in Fig. 3a.The time evolution in seismicity prior to a major earthquake is considered as a process of the earthquake fault system approaching the limit of spinodal stability.The characteristic earthquake eventually occurrs at the spinodal point (Rundle et al. 2000c).For an earthquake cycle with a span of T, the whole seismogenic process can be divided into three different periods, i.e. 1.0 > Dt > 0.2, 0.2 > Dt > 0.1, and 0.1 > Dt > 0.0.Here Dt = 1 -t T is the non-dimensional time remaining until the next characteristic earthquake and t is the time measured forward from the previous characteristic earthquake.During the entire cycle small earthquakes occur uniformly.By contrast, there is a systematic deficit in the number of moderate events during the earlier period 1.0 > Dt > 0.2, which is followed by an increase in the number of moderate events in the mid period 0.2 > Dt > 0.1.Seismic activation and further increase in the number of moderate events occurs in the precursory period 0.1 > Dt > 0.0.Frequency-magnitude distributions of earthquakes in the study area before the Pingtung doublets exactly mimic the seismicity evolution described above.Following the correlation-coefficient calculation of FMDs as proposed in Chen ( 2003), we obtained two characteristic times, i.e., March 2002 andMarch 2004, separating the earthquake cycle.We observed a deficit of moderate events with magnitude 5 ~6 before March 2002 (Fig. 3c) and a gradual increase in the number of moderate events from March 2002 to March 2004 (Fig. 3d).After March 2004 and before the Pingtung doublets the number of moderate events increased further (Fig. 3e), indicating the precursory activation of moderate earthquakes in the study area.Carefully inspecting the vertical intercepts at magnitude 3.2 in Figs.3c ~e, it seems there exists an increasing trend in the earthquake numbers for those three periods from Fig. 3c to e.We thus counted the monthly number of earthquakes occurring in the study area, and again found an increasing trend in the time series of monthly earthquake numbers (Fig. 3b).Two possible choices of t 1 have been determined from the FMD analysis above.The PI maps considering the result of the FMD analysis are shown in Fig. 4a   surrounding the hypocentral area of the Pingtung doublets, rather than a concentration at the hypocentral area.Two PI maps in Figs.4a and b may suggest regions in which intense seismicity changes gradually migrated from the surrounding area to the hypocentral area of the Pingtung doublets as the characteristic main shock approached.We will postpone the investigation of the migration pattern for the Pingtung doublets to a future work.

DISCUSSION
The tectonic setting probably affects the PI calculation.Different PI maps calculated from different depth selections are also shown in Figs.4a, c, and d, for which the ranges are 30 ~100 km, 0 ~100 km, and 0 ~30 km, respectively.In Fig. 4d it is quite reasonable that no sign for the Pingtung doublets could be found since events deeper than 30 km were not involved in the PI calculation.On the other hand, when comparing Fig. 4c with Fig. 4a, the hotspots around the Pingtung doublets become unclear in Fig. 4c.Clearly, the seismic activity is much higher within the depth from 0 to 30 km, as shown in Fig. 1c, and a strong signal arising from the shallower region can dominate the hotspot distribution.Therefore, the hotspots in Fig. 4c lie in almost the same locations as in Fig. 4d and the sign for the Pingtung doublets become contaminated and weak.These results suggest that mixing events of different depth scales and earthquake fault systems in different tectonic setting is problematic, and might have to be separately calculated in the context of the PI method.
Since the PI map is derived by computing the seismicity rate changes in the change interval relative to the background seismicity, the change interval plays an important role in the PI calculation.It is enlightening to compare this study with the PI analysis of the 1999 M L 7.3 Chi-Chi (Taiwan) earthquake (Chen et al. 2005).In this study, the earthquake fault system began to show major seismicity changes beginning in March 2004 (Fig. 1e).The change interval is about three years in the present case of the Pingtung doublets, but it is almost six years for the Chi-Chi earthquake (Chen et al. 2005).It seems reasonable that the preparation process is longer for a larger earthquake than for a smaller one.Such a speculation is also supported by the PI analysis, conducted by Chen et al. (2006b), for the 1992 M L 7.2 Landers (California) earthquake.The Landers and Chi-Chi earthquakes have similar magnitudes and rupture lengths of about 100 km and also, the preparation processes were almost six years in both cases.

CONCLUSION
Hotspots determined by the PI method indicate locations of anomalous seismic activity, associated with the occurrence of upcoming large earthquakes.The precursory period of anomalous seismicity associated with seismic activation can be determined by investigating the FMD based on the SOS model (Rundle et al. 2000c;Chen 2003).In the case of the Pingtung doublets, the precursory activity in the hypocentral area was seen to begin in March 2004 (Fig. 3e), with hotspots on the PI map clearly showing anomalous regions around the hypocenter (Fig. 4).In our approach, the FMD analysis aids the determination of the change interval in the PI calculation.By inspecting temporal variations in the FMD of earthquakes and the hotspot locations with intense seismicity changes on the PI map, it is possible to reveal the relevant preparation processes of the forthcoming catastrophic earthquakes.

Fig. 1 .
Fig. 1.(a) Epicenter distribution (blue circles) of seismicity used in this study.Two stars represent the M L 6.7 (21.6783°N, 120.5533°E) and M L 6.4 (21.9698°N, 120.4179°E)Pingtung earthquakes on 26 December 26 2006.(b) Distribution of cumulative number of earthquakes versus magnitude (blue line).The red line is the Gutenberg-Richter relation fitted from events with magnitude larger than 3.2.(c) Histogram of hypocentral depth for earthquakes.Bin width is 10 km.
can be computed and denoted by DI (x i , t b ). 4. The intensity change DI (x i , t b ) is an explicit function of x i and t b .For each x i , shifting t b produces a time series of [DI j (t b with fixed x i )].We then carry out a temporal normalization, i.e., subtracting the temporal mean [DI (t b with fixed x i )] and dividing by the temporal standard deviation s [DI (t b with fixed x i )], and get a temporally normalized intensity change D Ĩ (x i , t b ).Similarly, for each t b , we spatially normalize the space series of [D Ĩ j (x i with fixed t b )] and consequently get a spatiotemporally normalized intensity change D $ I (x i , t b ). 5. The normalized intensity change D $ I (x i , t b ) is computed at each location and its absolute value is further taken.The temporal average of DÎ (x i , t b ) , denoted by D $ ( )

Fig. 3 .
Fig. 3. (a) Frequency-magnitude distributions of three distinct periods for a characteristic earthquake cycle based on the self-organized spinodal model (Rundle et al. 2000c).For the detailed explanation please refer to Section 4. (b) Monthly number of earthquakes with magnitude larger than 3.2 occurred in the study area.Red dash lines mark March 2002 and March 2004, respectively.Black dash arrow indicates a vague trend of increasing event numbers.Also appearing are some peaks representing bursts of activity around June 2000 and March 2002.(c) ~(e) Normalized frequency-magnitude distributions in three periods as proposed in (a).There exists a systematic deficit in the number of moderate-size events during the earlier period (c), followed by an increase in the number of moderate events in the mid period (d).Seismic activation and further increase in the number of moderate events occurred in the precursory period (e) from March 2004 through the time of the Pingtung earthquake.Red lines in (c) ~(e) all have the same slope as the fitted Gutenberg-Richter relation in Fig. 1b.
Fig. 4. PI maps obtained from four choices of PI parameters: (a) t 1 = 2004/03/01 and hypocentral depth between 30 and 100 km; (b) t 1 = 2002/03/15 and hypocentral depth between 30 and 100 km; (c) t 1 = 2004/03/01 and hypocentral depth between 0 and 100 km; (d) t 1 = 2004/03801 and hypocentral depth between 0 and 30 km.Colored pixels (hotspots) represent areas with intense seismicity changes caused by both the seismic activation and/or quiescence.Blue circles represent earthquakes occurring in a short span of six months after t 2 (= 2006/05/01), with magnitude larger than 4. Note that the hypocentral depths of Pingtung earthquakes are outside the selected depth range in the PI map (d), thus two stars are added to indicate the epicenter locations of the Pingtung doublets.