Preliminary Results of the iSTEP Program on Integrated Search for Taiwan Earthquake Precursors

1 Institute of Geophysics, National Central University, Chung-Li, Taiwan, ROC 2 Institute of Space Science, National Central University, Chung-Li, Taiwan, ROC 3 Center for Space and Remote Sensing Research, National Central University, Chung-Li, Taiwan, ROC 4 Institute of Statistics, National Central University, Chung-Li, Taiwan, ROC *Corresponding author address: Prof. Yi-Ben Tsai, Institute of Geophysics, Chung-Li, Taiwan, ROC; E-mail: ybtsai@geps.gep.ncu.edu.tw The predictability of earthquakes has been a hotly debated question in earthquake science for some time. The answer to which begs another question, “Are there credible earthquake precursors?” Intrigued by these questions and encouraged by instrumental observations of conspicuous ionospheric and geomagnetic disturbances before the disastrous 1999 Chi-Chi earthquake, we decided to undertake an integrated Search for Taiwan Earthquake Precursors, called the iSTEP program. The multidisciplinary program includes five major components aimed at identifying potential seismological, geomagnetic, geodetic and ionospheric precursors, and statistical testing of any identified precursors. Since the program’s inception in April 2002, some encouraging preliminary results have been obtained. These includes precursory P wave travel-time changes six years before the Chi-Chi earthquake, identifiable geomagnetic changes two years before M ≥ 6 earthquakes, consistent ionospheric anomalies three days before M ≥ 5 earthquakes. We have also developed high-resolution radar interferometric methods for monitoring crustal deformation. A method for forecasting aftershock distribution on the basis of stress transfer has also been successfully tested on the Chi-Chi earthquake sequence. Aided by the newly installed geomagnetic and ionospheric networks we are hopeful about identifying more earthquake precursors and developing relevant physical mechanisms in the remaining two years of the iSTEP program. (


INTRODUCTION
In a region with high seismicity and dense population, Taiwan has suffered numerous disastrous earthquakes in the past.The steady collision between the Philippine Sea and Eurasian plates means that Taiwan will inevitably face earthquake hazards in the future.The 1999 Chi-Chi earthquake, which took more than 2,500 lives and destroyed more than 100,000 dwellings, was only the latest case.For some natural hazards, such as typhoons, short-term forecasting is an integral component of preparedness.In response, buildings can be secured, equipment can be removed, emergency services can be put on alert, and populations can be evacuated, if necessary.Unfortunately, earthquakes still cannot be forecast, much of this is due to a lack of predictive capability.
The problem of earthquake prediction has been under intensive investigation for nearly forty years in Japan, China, and the U.S. Yet relatively little progress has been made thus far.Long frustration over the surprising difficult area of earthquake prediction has recently led to serious debate within the seismological community regarding the predictability of earthquakes.Some people argue that earthquakes are inherently unpredictable (e.g., Geller 1997).Yet others believe that earthquake prediction is possible (e.g., Whiteside 1998).Our present inability to predict earthquakes is partly due to inherent characteristics of earthquakes and seismic waves and partly to incomplete understanding of the earthquake process (Silver 1998).At the heart of the debate is a question of whether or not there are recognizable precursors of earthquake.
Many precursors have been mentioned, but it is not clear which, if any, are reliable.Certainly any reliable scheme for practical prediction must be based on a combination of clues, so that firm decisions are possible before warnings are issued.Several promising clues have been mentioned in the literature, such as detection of strain in the Earth's crust by geodetic surveys, the identification of suspicious gaps in the regular occurrence of earthquakes in both time and space, and the observation of foreshocks.In recent years, major earthquake prediction research has emphasized more precise measurements of fluctuations in physical parameters of crustal rocks in seismic areas.According to the dilatancy-diffusion model (Scholz et al. 1973), the five parameters thought to be relevant are: seismic P velocities, the uplift and tilt of the ground, the emission of radon gas from water wells, electrical resistivity in rocks, and the number of earthquakes in the region.
An earthquake prediction research program in Taiwan was initiated in 1977, jointly by Y. B. Tsai of the Academia Sinica and T. L. Teng of the University of Southern California, and funded separately by the National Science Council and the U.S. Geological Survey (Tsai et al. 1983).Under this earlier program all the above five parameters were observed.Some of the parameter observations, such as radon emission were discontinued due to funding and manpower constraints.Others were upgraded to more advanced methods, such as replacement of traditional trilateration measurements of crustal deformation with the new global positioning system (GPS) methods.Still there are other areas that have continued using original instruments until this day for observing changes in geomagnetic intensity and seismicity.
Recently, three new developments have prompted us to undertake a new program to conduct integrated search for Taiwan earthquake precursors, called iSTEP.Fortuitously, clear changes in geomagnetic intensity were observed at one of the original stations at Liyutan at least one month before and after the 1999 Chi-Chi earthquake (Yen et al. 2004).This was the first major new development which encouraged us to undertake the present iSTEP program.The second major new development prompting us to undertake the present program was the observation of clear electromagnetic precursors in the ionosphere, both by ionosondes and GPS data several days before large earthquakes (Liu et al. 2000(Liu et al. , 2001)).
The third major new development was the physical process proposed recently (Freund et al. 1999;Freund 2000;Freund et al. 2004) to explain electromagnetic signals associated with earthquakes observed both on the ground and in the ionosphere.According to this process, igneous rocks contain positive hole pairs (PHP), i.e., dormant electronic charge carriers.The PHP can be activated by microfracturing and/or dislocation movement.Upon activation, the PHP release highly mobile charge carriers in the form of positive holes.Positive holes are defect electrons in the oxygen sublattice of minerals that can conduct through rocks.Once generated these charges spread out from the source volume.Further understanding of these positive hole charge carriers and their manifestation will enable us to reevaluate electrical phenomena associated with earthquake activity.
Considering the above background, the iSTEP program is set to achieve the following three objectives: 1.To identify any possible, yet reliable precursory phenomena of earthquakes, including seismological, geodetic, geomagnetic, and ionospheric variations.
2. To characterize feasible mechanisms behind the identified precursory phenomena in order to give quantitative physical interpretations.
3. To correlate the identified precursory phenomena with existing or emerging earthquake preparation models in order to lay a foundation for future earthquake prediction.

MAJOR COMPONENTS OF THE iSTEP PROGRAM
Given the above scientific objectives and taking into account the strength of our present research team, we decided to include the following five major components in the iSTEP program (Fig. 1).

A. Potential Earthquake Precursors Based on Seismological Variations
We use the excellent earthquake catalogs and high-quality digital seismic waveform data obtained by the Central Weather Bureau (CWB) to study the following topics.
1. Precursory changes in P wave velocities and waveforms.If rock properties change before an earthquake, then the P wave velocity and waveforms may also vary.We will assemble the P wave residuals and spectral contents at each of the Central Weather Bureau Seismic Network (CWBSN) stations as function of time to see whether there are clear changes (i.e., delays) preceding large earthquakes of M ≥ 6.0 in the nearby areas.
2. Variations in seismicity rate: Significant changes in normal background seismicity were observed for several large earthquakes.Another method is to measure the changes in b value.We will perform a complete search of all M ≥ 6.0 earthquakes to identify possible seismicity rate changes.
3. Changes in shear wave splitting: Shear wave splitting can be used to monitor the effects of stress build-up before earthquakes, and hopefully to forecast future large earthquakes (Zatsepin and Crampin 1997;Crampin and Zatsepin 1997;Crampin 1998).The goal of this study is to systematically examine possible anisotropy changes before earthquakes and look for possible earthquake precursors.Especially, we explore possible correlation of these changes with the electromagnetic disturbances in ionosphere.
4. Stress transfer among large earthquakes and the physics on earthquake dynamic triggering will be addressed.The correlation among these observations will be examined.Seismological observations will further link with other non-seismological observations from other iSTEP program components to find possible earthquake precursors and to look for feasible physical mechanisms of these observations.

B. Potential Earthquake Precursors Based on Variations of Geomagnetic and Gravity Fields
After the Chi-Chi earthquake, we examined data from a geomagnetic network of eight stations operational since early 1980's to provide continuous observations of total intensity over Taiwan.It was found that variations in the eight stations generally yield a similar tendency.At the Lunping station, which served as the base-station of the geomagnetic network, no disturbance was recorded before and after the 1999 Chi-Chi earthquake.However, during the two months prior to the earthquake, data recorded at Liyutan geomagnetic station (24.2°N, 120.4°E), which is only about six kilometers north of the northern end of the Chelungpu fault, fluctuated significantly.It was surprising to find that these fluctuations disappeared on 22 October 1999, when a strong earthquake M = 6.2 occurred near the southern end of this fault.These anomalous variations, more than 150 gammas in amplitude, appeared to be associated with the earthquake (Yen et al. 2004).As far as we know, this type of geomagnetic change has never before been reported in scientific literature.Accordingly, we include geomagnetic variations as part of the iSTEP program.
We will first analyze the continuous geomagnetic data observed over the last 20 years at the eight permanent stations.We will also upgrade the existing geomagnetic network (Fig. 2).In addition to geomagnetic methods, repeated microgravity surveys along profiles across known active faults will also be made to detect crustal deformation, subsurface material movement and subsequent mass redistribution.

Interferometry
The use of interferometric synthetic aperture radar (InSAR) for monitoring crustal deformation processes has received considerable attention lately.Several studies suggested great potential for the InSAR technique for mapping coseismic movement associated with strong earthquakes (Massonnet et al. 1993;Zebker et al. 1994).Microwave frequency allows for surface changes to be detected on a scale of centimeters.A number of studies have demonstrated the potential and capability of satellite radar interferometry for mapping ground surface deformation.For small-scale surface height movements, differential radar interferometry is usually adopted.
The objectives of this iSTEP program component can be summed up as follows: 1.To integrate atmospheric and ionospheric observations from very high frequency (VHF) and GPS into radar interferometry for path length calculations, and thus achieving higher estimation accuracy of terrain height and height changes.
2. To detect surface movements and pre-, co-, and post-seismic crustal deformations.3. To analyze and investigate the correlation of radar interferograms with the temporal variations and spatial distributions of seismic electromagnetic phenomena in the ionosphere and atmosphere.
4. To process and analyze the thermal infra-red (IR) satellite imagery as a precursor before an earthquake.

D. Potential Earthquake Precursors Based on Ionospheric Variations
Although seismic waves are the most obvious manifestation, earthquakes may be accompanied or preceded by signals of a different nature -electric, electromagnetic, or luminous -that may help forecasting impending seismic activity (Uyeda 2004;Freund 2000).After the disastrous 1995 Kobe earthquake in Japan, the Science and Technology Agency of Japan called for the "Earthquake Integrated Frontier Research" project.It started in 1996 to study extensively the whole spectrum of seismo-phenomena taking place in the atmosphere and ionosphere/magnetosphere and then to investigate energy transfer from the lithosphere to the atmosphere and ionosphere/magnetosphere (Uyeda 2004).We analyzed ionospheric foF2 and found three clear precursors appearing about 1, 3, and 4 days prior to the Chi-Chi earthquake.In this iSTEP program component, we will systematically search through ionospheric variations of all large earthquakes in Taiwan to see whether there are persistent precursors.Both ionosonde and GPS data will be used (Fig. 2).

E. Statistical Studies of Potential Earthquake Precursors
Although extensive and reliable data on possible precursors of major earthquakes are expensive to collect, and may take a long time to accumulate (Vere-Jones 1995), we have observed some surprising anomalous fluctuations of the maximum plasma frequency and total electron content (TEC) of the ionosphere which may possibly be related to strong earthquakes in Taiwan area (Liu et al. 2000(Liu et al. , 2001(Liu et al. , 2004a)).These anomalies reveal not only the temporal, but also the spatial information about the earthquakes to follow.Therefore, according to the empirical data we have obtained so far, we will develop appropriate statistical process control (SPC) techniques (Stoumbos et al. 2000) to identify the variations of the relevant measurements in electromagnetic fields.
In order to validate the reliability of precursory parameters, the Monte Carlo simulations (Musson 1997;Stark 1997) of some reasonable prediction processes based on either the occurrence of precursors or previous earthquake catalogs will be compared with the earthquakes of interest.The generalized linear model (McCullogh 2000) is further introduced to describe the relationship between the occurrence or strength of precursors and the associated earthquake characteristics, for instance, magnitude, depth and the distance of propagation from the epicenter to the ionosonde station.

PRELIMINARY RESULTS OF iSTEP PROGRAM
The iSTEP program started officially in April 2002 for four years.Since then an upgraded geomagnetometer network has been installed.Several ionosonde stations have also been added.Figure 2 shows the locations of these instruments.Since data obtained by the new observational facilities are still in the early stage, we concentrated our research efforts in the first two years on analyzing the existing data for past earthquakes.Detailed results are presented in several papers both in this special issue and elsewhere (Lee and Tsai 2004;Chan and Ma 2004;Chen et al. 2004;Chang et al. 2004;Liu et al. 2004b;Chen et al. 2004).We will only give a summary below to highlight the results obtained so far.

A. Seismological Results
1. P wave velocity varitations before and after the Chi-Chi earthquake P wave velocity variations may provide precursory information before large earthquakes.We present results on variations of P wave travel-time residuals before and after the 1999 Chi-Chi earthquake (Lee and Tsai 2004).The short-period seismic network deployed by the CWB (CWBSN) is routinely used for earthquake location in Taiwan.Since the CWB used a 1D velocity model, the travel-time residuals would show lateral and vertical inhomogeneities of the velocity structures.Also, the travel-time residuals may change with time if the velocity structure changes temporally.We used the CWBSN data from 1991 to 2002 to obtain the P wave travel-time residuals.Our purpose is to find the mean value of the P wave travel-time residuals at each station and its temporal changes at stations near the Chelungpu fault.
The results show that the mean P wave travel-time residuals at stations 40 km east of the Chelungpu fault changed very little both before and after the Chi-Chi earthquake.But the mean P wave travel-time residuals at the stations immediately west of the Chelungpu fault changed significantly before and after the Chi-Chi earthquake.The anomalous zone is bounded by stations in the black frame, as shown in Fig. 3.These temporal changes of P wave traveltime residuals can be attributed to changes in the velocity structure, which in turn might be caused by crustal deformation east of the Chelungpu fault beginning about six years before the Chi-Chi earthquake.In other words, there was a long-term six-year precursor of P wave velocity decrease before the Chi-Chi earthquake which may be associated with dilatancy due to development of cracks.

Forecasting aftershock distribution from stress changes following large earthquakes
Non-uniform spatial slip dislocation models from seismic waveform inversions of several past earthquakes in Taiwan are used to calculate possible stress transfer associated with the aftershock distribution.Our results show that the positive stress changes calculated for optimal orientation of fault planes after the mainshock are correlated with the aftershock distribution.In order to assess the possibility for forecasting aftershock distribution from stress changes of the mainshock, we used a homogeneous fault model on the basis of the earthquake scaling law to make rapid stress change calculations.
Comparison of the stress changes from both the homogeneous and heterogeneous fault models show similar patterns with good correlation with aftershock distribution, including the complex fault rupture of the Chi-Chi earthquake (Fig. 4).Our results, thus, indicate good possibility for forecasting aftershock distribution from stress changes of the mainshock.Once the location, magnitude and focal mechanism of a strong earthquake become available, stress change calculations can be carried out to forecast aftershock distribution for timely aftershock hazard mitigation.

B. Geomagnetic Total Field Variations Associated with Earthquakes in Taiwan
Electromagnetic phenomena associated with seismic activity have been extensively discussed (Pulinets 1998;Molchanov and Hayakawa 1998;Hayakawa 1999;Freund 2000).Zeng et al. (2001) studied a large numbers of cases and found about 80% of M ≥ 6.0 earthquakes occurred within nine months to 2.5 years after appearances of the zero isoporic zone (ZIZ), which is defined as the annual change rate of geomagnetic parameters to be less than or equal to ± 5 nT yr -1 .
We have examined variations in the geomagnetic total field recorded by the eight magnetometer stations which might be related to occurrence of M ≥ 6.0 earthquakes in Taiwan from 1989 to 2001.It is found that the annual change rates of the field strength reduced down to a zero isoporic value of ± 5 nT yr -1 after 1997 and the following M ≥ 6.0 earthquakes often occur within the ZIZ.Our results, as shown in Fig. 5, suggest that the spatial and temporal ZI signatures tend to lead the occurrence of M ≥ 6.0 earthquakes by about two years.

Deformation in Taiwan
The orogen of Taiwan is young and active as revealed by frequent earthquakes.Some strong earthquakes, with the accompanying crustal deformation often caused severe damages.We have reviewed recent results from application of the radar interferometry technique to monitor crustal deformation in Taiwan.Results from five deformation events have been obtained.They are the coseismic deformation of the Chi-Chi earthquake, uplift of the Tainan area, active deformation of the Hukuo area, rapid land subsidence of the Chungli area, and seasonal land subsidence in the Pingtung plain's area.As shown in Fig. 6, for the uplift in the Tainan area, the results show that the radar interferometry is a useful high-resolution tool for monitoring crustal deformation of different characteristics.

D. Ionospheric foF2 and TEC Anomalies During M ≥ 5.0 Earthquakes in Taiwan
We have examined variations in the critical frequency foF2 recorded by an ionosonde, and the TEC derived from a network of 5 ground-based receivers of the GPS, and correlated them with the occurrences of 144 M ≥ 5.0 earthquakes in Taiwan during 1997 -1999.Results in Fig. 7 show that the foF2 and TEC yield similar patterns, and often concurrently register pronounced depression anomalies four days before the earthquakes.An unbiased investigation of all anomalies before and after the earthquakes positively confirms that the anomalous depression in the foF2 and TEC are mostly premonitory anomalies.

E. Statistical Tests on Pre-earthquake Ionospheric Anomalies
Anomalous depression of the maximum plasma frequency in the ionosphere (foF2) appears consistently within 1 -5 days before many M ≥ 5.0 earthquakes in the Taiwan area during 1994 -1999 (Liu et al. 2003), as shown in Fig. 8.We have made two statistical tests

CONCLUSIONS
The iSTEP program is aimed at searching for earthquake precursors using multidisciplinary approaches.Since its inception in April 2002, some preliminary results have been obtained and summarized below.
From seismological studies, we obtained two main results.First, we examined the P wave velocity variations before and after the Chi-Chi earthquake by analyzing the P wave traveltime residuals.The result shows P wave velocity changes in the footwall of the Chelungpu fault beginning about six years before the Chi-Chi earthquake.This appears to signal precursory phenomena due to development of cracks and resultant dilatancy of rocks.Second, a model for forecasting aftershock distribution from stress transfer following large earthquakes was tested by the Chi-Chi earthquake.The results indicate good possibility for forecasting aftershock distribution from stress changes due to the mainshock.If we have sufficient timely information of a large earthquake, we can forecast its aftershock distribution for effective aftershock hazard mitigation.
In geomagnetic studies, we examined variations in the geomagnetic total field recorded by eight magnetometer stations which may be related to occurrence of M ≥ 6.0 earthquakes in Taiwan.The results show that the annual change rates of the field strength reduced down to a zero isoporic value of ± 5 nT yr -1 after 1997.The following M ≥ 6.0 earthquakes often occur within the ZIZ.Moreover, the spatial and temporal ZI signatures tend to lead the occurrence of M ≥ 6.0 earthquakes by about two years.
By processing the radar interferometry images, we observed several types of ground deformation, i.e., uplift of the Tainan area, active deformation in the Hukuo area, rapid land subsidence of the Chungli area, and seasonal land subsidence of the Pingtung plain's area.Thus radar interferometry is demonstrated to be a useful high-resolution tool for monitoring crustal deformation of different characteristics, including that of tectonic origin.
From ionospheric observations, we examined variations in terms of foF2 and TEC and correlated these variations with the occurrence of earthquakes with M ≥ 5.0.The variations of both foF2 and TEC show similar patterns, and often concurrently register pronounced depression anomalies four days before large earthquakes.Unbiased statistical tests of all anomalies before and after the earthquakes positively confirm that the anomalous depression in the foF2 and TEC are mostly premonitory anomalies.
Two statistical tests were made against the foF2 anomaly as a candidate precursor of earthquakes.One is designed to compare the foF2 anomaly-based method with competitive alternatives for predicting the earthquakes under study.The other is to investigate the significance of the identified foF2 anomalies, among all possible foF2 anomalies, related to the recorded earthquakes with M ≥ 5.0.The results show that the foF2 anomaly is superior to the alternatives for temporal prediction of the M ≥ 5.0 earthquakes in Taiwan area.
The results summarized above, even though still preliminary, appear to have identified several promising precursory changes, ranging from a few years to a few days prior to large earthquakes.In the case of forecasting aftershock distribution, we in fact treat a mainshock as a certain precursor.The iSTEP program will continue for at least two more years.We hope to identify more earthquake precursors from observations of our new instrumental networks in the near future.In the meantime, we will pursue relevant physical mechanisms.Similar programs are in progress in other countries, such as Japan, China, and the U.S. International cooperation with these and other countries are important to us.Through integrated search for earthquake precursors and positive establishment of related physical mechanisms, we hope to lay the foundation for eventual earthquake prediction.

Fig. 1 .
Fig. 1.The five major components of the iSTEP program.

Fig. 2 .
Fig. 2. Locations of new magnetometer and ionosonde stations deployed under the iSTEP program.

Fig. 3 .
Fig. 3.A contour map of mean P wave travel-time residuals in Taiwan for three different time periods: (a) Mean P wave travel-time residuals from 1991 to 1993.(b) Mean P wave residuals from 1994 to Chi-Chi earthquake.(c) Mean P wave residuals from Chi-Chi earthquake to 2002.(After Lee and Tsai 2004).

Fig. 4 .
Fig. 4. The changes of Coulomb failure stress based on (a) the heterogeneous slip model by Ji et al. (2003), (b) three-segments model, and (c) onesegment model, respectively.The three-month aftershocks (green circles) are considered.The amount of Coulomb stress change is shown by the colored bar.(After Chan and Ma 2004).

Fig. 5 .
Fig. 5. Temporal variations of the geomagnetic total field during 1988-2001.The solid and dashed lines denote the yearly value of the observed and IGRF field, respectively.The error bars are the standard deviation of the yearly values.(After Chen et al. 2004).

Fig. 8 .
Fig. 8. Cross-correlation coefficient of the foF2 anomalies as a function of time relative to the origin time of earthquakes.(After Liu et al. 2003).