REVIEW A review on precursors of the 1999 M w 7.6 Chi-Chi, Taiwan, earthquake

Earthquake prediction has been a long-term debatable problem in earthquake science. There were numerous studies of possible precursors after the M w 7.6 Chi-Chi, Taiwan, earthquake of 20 September 1999. A deep study about these precursors will be useful for the exploration of the debatable problem. In this study, based on the time window (or the precursor time, T ) prediction for earthquakes in Taiwan are defined as follows: very-long-term prediction ( T > ten years or longer); long-term prediction (ten years> T > three years); intermediate-term prediction ( T = six months to three years); short-term prediction ( T = eight days to six months); and imminent prediction ( T ≤ seven days). Meanwhile, the precursors are classified into four cat - egories: mechanical, electromagnetic, geochemical, and biological precursors. Each category may include several items. First, we will compile all given possible precursors and the respective precursor times. Secondly, we will examine these possible precursors based on known physical and chemical theories. Thirdly, we will discuss the possible correlations between two precursors. From the results of this study, most of precursors are reliable and few of them need re-examinations. This suggests the possibility of earthquake prediction. Nevertheless, earthquake scientists cannot yet preciously predict an earthquake just based on the reliable observations. Like weather forecasting, it is very necessary to construct the physical and chemical models of respective precursors or even a unified model for all precursors. Based on the models or a unified model, it is possible to predict earthquakes.


INTRODUCTION
In order to reduce seismic hazards, one of the significant ways is to warm the occurrence of an impending earthquake to the government and the public from reliable earthquake precursors. Milne (1880) who was a professor of Tokyo Imperial University, Japan wrote: 'Even since seismology has been studied one of the chief aims of its students has been to discover some means which enable them to foretell the coming of an earthquake….' The earthquake rupture processes are considered to be preceded by complex precursors that might happen during the nucleation processes prior to an impending earthquake. When reliable precursor may be preciously detected, it is possible to give effective warning of the event. Previous studies (Turcotte 1991;Keilis-Borok 2002;Uyeda et al. 2009) show four different categories of precursors: mechanical, electromagnetic, geochemical, and biological precursors. For each category, there are several items. Recently, a book (Ouzounov et al. 2018) that was published by the American Geophysical Union, USA describes multidisciplinary approaches to earthquake prediction.
After the 1891 Nobi earthquake, Japanese founded the Imperial Earthquake Investigation Committee (Imamura 1937). Tsuboi et al. (1962) proposed the 'Blueprint' that was the basis for Japan's prediction program and is almost the first national program for earthquake prediction around the world. Long-term, extensive searching for reliable earthquake precursors done in the last one hundred years does not seem successful (Geller 1997;Sornette 1999; and cited references therein). Historically, there was a famous debate between Prof. Omori and Assit. Prof. Immamura of Tokyo Imperial University about the possible occurrence of a big earthquake in the Kanto area before 1923. The former considered that earthquakes are random and un-predictable; while the latter assumed that earthquakes are predictable. However, a M 7.9 earthquake, which caused serious damages, happened on 1 September 1923 in the Kanto area. Of course, the latter did not predicted the occurrence time and size of that event. In 1999, there was a forum for debating the problems of earthquake prediction was chaired by Prof. Main of Edinburgh University, England, and the debates were published in the journal 'Nature.' Clearly, earthquake prediction has been a paradox of earthquake science (Wang 2019) and the most difficult challenge of earthquake scientists. Aki (1989) essentially assumed that earthquakes are predictable. He formulated ideal prediction as follows: we would inform the public that the probability of the occurrence of an earthquake is one per T units of time (for example, years) before a target event with given magnitude and location. He also defined the task of a geoscientist in earthquake prediction as objectively estimating the probability of occurrence of an earthquake with a specified magnitude, place, and time window under the condition that a particular set of precursory data was observed. Aki (1995) emphasized the societal implications of earthquake prediction. Numerous geoscientists claimed the possibility of earthquake prediction (e.g., Rikitake 1975Rikitake , 1979Rikitake , 1984Rikitake , 1987Park et al. 1993;Knopoff 1996a, c;Wyss et al. 1997;Uyeda et al. 2009;Freund 2013;Varotsos et al. 2017). But, some others (e.g., Geller 1996Geller , 1997Geller et al. 1997;Kagan 1997;Leary 1997; and cited references therein) strongly argued the possibility of earthquake prediction. Geller (1997) described in details the history of debates on earthquake prediction. Geller (1997) and Geller et al. (1997) explained their reasons why earthquakes cannot be predicted.
There are two basic problems from the previous discussion. The first problem is the reliability of observed precursors, because numerous precursors have not been tested physically and statistically. The second problem is how to use the reliable precursors to predict an earthquake. For the first problem, several authors (e.g., Molchan and Kagan 1992;Kagan and Jackson 1996;Kagan 1997) stressed the importance of statistical tests on observations before using them to predict earthquakes. This means that the observed precursors must be tested through statistical methods (Kagan 1997). For instance, Kagan and Jackson (1996) tested fourteen earthquake predictions for M ≥ 5.8 events during 1987 to 1995 in Greece using VAN's methodology (Varotsos et al. 1981). They concluded that only one case for a swarm with the maximum magnitude of 5.8 passed the test, while the rest 13 events all failed. This is a very serious problem for earthquake prediction (Geller 1996(Geller , 1997. Molchan and Kagan (1992) suggested the idea and methodology for statistical tests of observed precursors.
The second problem is due to a lack of acceptable working models for respective precursors or a unified mod-el for all precursors. We cannot predict an impending earthquake just from observed data, even all observed precursors are reliable. Like weather forecasting, earthquake scientists have first to construct an acceptable working model based on reliable precursors and then use the model to predict earthquakes. Up to date, this is still a difficult challenge for earthquake scientists, even though there have been some physical models (cf. Wang 2016Wang , 2017Wang , 2018bWang , 2020Wang , 2021 and cited references therein).
For the precursor time, T, Wallace et al. (1984) and Kisslinger (1989) first used 'long-term' (T = a few years to a few decades), 'intermediate-term' (T = a few weeks to a few years), and 'short-term' (T < a few weeks) to show the windows of prediction. Of course, it is difficult to define the exact time windows that may be distinct in different seismogenic zones. Wang et al. (2018) defined four time windows for earthquake prediction in China: ultralong and long term time period, medium term time period, short term time period, and impending period. However, they did not clearly define the length of each time window. Here, I define five types of earthquake prediction, especially for Taiwan: verylong-term prediction, long-term prediction, intermediateterm prediction, short-term prediction, and imminent prediction. The time windows of the five types of prediction will be explained below. In fact, the very-long-term prediction is based on earthquake recurrence that inferred from historical documents and data from geological trenching surveys. Since the recurrence time is usually longer than several hundred years for large earthquakes, earthquake recurrence is more appropriate for forecasting than for prediction.
On 20 September 1999, the M s 7.6 Chi-Chi earthquake ruptured the Chelungpu fault, which is a ~100-km-long and east-dipping thrust fault, with a dip angle of ~30°, in central Taiwan (Ma et al. 1999;Shin 2000;Shin and Teng 2001). The epicenter (denoted by a star) and the surface trace of the fault (displayed by a solid line with 'CLPF') are shown in Fig. 1. (Included also in Fig. 1 are two nearby faults for references.) Unfortunately, this event was not predicted or forecasted by Taiwan's earthquake scientists because it occurred along the Chelungpu fault that was classified to be the Type-II active fault before 1999. After the earthquake, numerous earthquake scientists examined the recorded data and tried to find out the possible precursors. All results used in this study were indeed made after the earthquake. Tsai et al. (2006Tsai et al. ( , 2018 summarized some precursors that are denoted by a symbol ' * ' in Table 1 and also schematically showed the precursor times. However, a few precursors shown in their study are questionable and the results done by some earthquake scientists were not included in their papers. In this study, a complete set of possible precursors for the earthquake is compiled and discussed in details. The very-long-term prediction is related to earthquake recurrence inferred from geological data and historical documents (Chen et al. 2001a(Chen et al. , b, 2004bOta et al. 2001; Wang 2005). The long-term prediction include the mechanical precursors: temporal variation in seismicity patterns (Chen 2003;Chen et al. 2005aChen and Wu 2006;Wu and Chiao 2006;Wu and Chen 2007;Wu et al. 2012), b-value anomalies (Tsai et al. 2006), and changes in P-wave traveltime residuals (Lee and Tsai 2004); the intermediate-term prediction include the mechanical precursors: pre-seismic crustal deformations (Yu et al. 2001), seismic quiescence (Wu and Chiao 2006;Wu and Chen 2007), and groundwater level changes (Chen et al. 2013(Chen et al. , 2015. The short-term prediction include (1) mechanical precursors: slow slip (Lin 2012); (2) electromagnetic (EM) precursors: geomagnetic anomalies (Yen et al. 2004;) and ULF signals (Akinaga et al. 2001); (3) geochemical precursors (Song et al. 2003;Wang et al. 2005); and (4) biological anomalies for some animals . The imminent prediction includes (1) mechanical precursors: infrasound (Xia et al. 2011); (2) EM precursors: the anomalies of ionospheric total electron content (TEC) , f o F 2 anomalies (Chuo et al. 2002), ELF/ULF signals (Ohta et al. 2001), lightning (Tsai et al. 2006), and earthquake lights ; and (3) biological anomalies for some animals .
Through interview with local people living in the source area of the Chi-Chi earthquake, Chen et al. (2000) compiled numerous records about some possible precursors, including pre-seismic phenomena (anomalous animal activities, earthquake light, changes of wind, and abnormal skylight) and co-seismic phenomena (earthquake light, ab-normal sounds, abnormal gaseous odors, and initial motions). The anomalous animal activities could be traced back almost two months before the mainshock. The frequency of anomalous animal activities was high and sharply increased on the day just before the earthquake.
The studies of precursors of the Chi-Chi earthquake seem able to give us an opportunity to explore the possibility of earthquake prediction. Three working steps that will be made in this study are described as follows: First, we compile all given possible precursors and their precursor times. Secondly, we examine these precursors based on known physical and chemical theories. Thirdly, we explore the possible correlations between two precursors or among several precursors. From the results, we may tackle the debatable problem if earthquake prediction is possible or not.

TWO BASIC MODELS
In order to examine the observed precursors of the 1999 Chi-Chi earthquake, we need two basic physical models of earthquakes that are described below.

Reid Elastic Rebound Model
Reid's elastic rebound model (Reid 1910) is shown in Fig. 2. The width of deformed zone, 2w, within which the shear strain develops progressively across the fault prior to the earthquake, is μD/Δσ where μ is the rigidity of fault-zone material (Turcotte and Schubert 1982). Since the

Types of Prediction Precursors T Remarks
Long-term

Mechanical Precursors
Surface deformations* 3 years Tsai et al. (2006) Appearance of seismic quiescence 9 months Wu and Chiao (2006); Wu and Chen (2007) Change of groundwater level 250 days Chen et al. (2013Chen et al. ( , 2015 EM Precursors  Liu et al. (2001Liu et al. ( , 2004a; Chuo et al. (2002) Lightning* 4 days Tsai et al. (2006) ULF/ELF signals 4 hours Ohta et al. (2001) Earthquake light few hours Chen et al. (2000) Geochemical precursor 2 days Song et al. (2006) Animal anomalies (for some animals) ≤ 7 days Chen et al. (2000)  Note: *: The symbol denotes that the related article was cited by Tsai et al. (2006). values of Δσ for large earthquake are in general 1 -10 MPa, those of 2w are 5 -50 km. Since the model was originally based on the San Andrea fault, which is a nearly vertical, strike-slip fault, it cannot completely describe the faulting properties (mainly thrust faulting) of the Chelungpu fault. Nevertheless, it can still give a first-order approximation to the behavior of the fault, because it had a strong strike-slip component, especially on its northern segment. For the Chelungpu fault, the average displacement and the average static stress evaluated by several authors are, respectively, D = 3.1 -6.0 m and Δσ = 4.2 -10.0 MPa (Wang 2006). In order to estimate the value of 2w, we consider two values of μ. The first one is μ = 19 GPa measured at a borehole having a depth of 1100 m ). This leads to 2w = 5.9 -27.1 km, with an average of 16.5 km. Hence, the average width of deformed zone on one side of the fault is 8.3 km. The second one is μ = 30 GPa for common crustal rocks (Turcotte and Schubert 1982). This results in 2w = 9.3 -42.8 km, with an average of 26.1 km. Hence, the average width of deformed zone on one side of the fault is 13.0 km. Since the width and dip angle of the ruptured fault plane are, respectively, 20 km and 30° (Wang 2006), the length of projection of the fault plane on the ground surface is 17.3 km. From the spatial distribution of co-seismic ground surface displacements measured by Yu et al. (2001), large displacements were observed within 20 km to the east of the fault. Hence, the two kinds of evidence suggest that the estimate of the width of deformed zone seems more acceptable from μ = 30 GPa than from μ = 19 GPa, because the bottom of the fault is much deeper than the borehole and thus the rigidity of fault-zone rocks should be stronger at depths than at the shallow part.

Temporal Variations in Stress and Slip on a Fault
Before an Earthquake The temporal variations in stress σ and slip u on a fault (Atkinson 1984;Rudnicki 1988;Main and Meredith 1989) may be schematically displayed in Fig. 3 (a thick solid line for σ and a thin solid line for u). A detailed description about the temporal variations in σ and u can be seen in Main and Meredith (1989) and Wang (2021). Figure 3 shows that the stress increases with time and then decreases with increasing time after it reaches its maximum value at t = t m ; while the slip increased very slowly with time from t 0 to t c before step 3c and then increased rapidly in step 3c up to the failure point at t r . The stress and slip are two important factors in influencing the generation of precursors (e.g., Dieterich 1978Dieterich , 1979Main 1988;Scholz 1990).
Based on the stress, there are three stages, i.e., Stages 1 -3, before an earthquake (at t = t r ) and two stages, i.e., Stags 4 -5, during and after the event. In general, the time period is the longest (several ten or hundred years) in Stage 1 from t o to t y and the shortest (usually several ten seconds) in Stage 4. In Stage 1 elastic buildup of strain energy starts at t o and stops at t y when the elastic (linear) loading transfers to plastic (nonlinear) loading, thus yielding the long-term precursors such as changes of seismicity patterns and the b-value anomalies. The changing point from elastic loading to plastic loading is called the yielding point. Stage 2 (occurring from t y to t m when the loading stress reaches its peak) displays plastic strain hardening due to dilatancy because of crack coalescence and fluid transport in the fault zone. Intermediate-term precursors, for examples, groundwater level change and hydrochemical anomalies, may be observed from this stage. In Stage 3 (occurring from t m to t r ) the loading stress decreases and demonstrates precursory stress drop or strain softening. Stage 3 has three steps (3a, 3b, and 3c): (3a) for microcrack linkage, (3b) for pore fluid diffusion, and (3c) for quasi-static slip on the fault between asperities. The mechanic behavior in this stage could result in short-term and imminent precursors such as the non-volcano tremors (e.g., Katakami et al. 2018), foreshock activity, infrasound, and EM signals. Figure 3 shows that the slip u increased very slowly with time before step 3c and then rapidly increases with time within this step until the failure of an event. In Stage 4 the fault breaks and an earthquake occurs at time t r , thus resulting in dynamic slip of the fault behind the crack tip. In Stage 5, aftershocks may last for several months or even several years. The possible time windows of long-term (starting from t 0 ), intermediate-term (starting from t y ), short-term (starting from t m ), and imminent (starting from t c ) precursors are also displayed in Fig. 3. Of course, the time window of very-long-term prediction is much longer than t 0 -t r and thus not shown in Fig. 3.
The weather bureau or agency of numerous nations usually provides weather forecasting in coming seven days to the public. Hence, based on weather forecasting and the occurrence times of precursors of the Chi-Chi earthquake, I take a time interval of 'seven days' to be an optimum time window for imminent prediction. Hence, I define five categories of prediction having different time windows for earthquakes in Taiwan as follows: very-long-term prediction (T > ten years or longer); long-term prediction (ten years > T > three years); intermediate-term prediction (T = six months to three years); short-term prediction (T = eight days to six months); and imminent prediction (T ≤ seven days).

OBSERVATIONS
Except for the earthquake recurrence, the values of T for respective precursors are listed in Table 1. Since the anomalous animal activities  occurred in both the short-term and imminent time windows, related observations will be placed on an independent section.

Very-Long-Term Prediction: Earthquake Recurrence
The recurrence time, T R , of a fault increases with earthquake magnitude, M, and its value may be several hundred or even a few thousand years. Hence, it is necessary to use historical documents and paleoseismicity inferred from geological trenching data to estimate the recurrence time. Numerous trenching surveys were conducted on the Chelungpu fault (Chen et al. 2001a(Chen et al. , b, 2004b(Chen et al. , 2005b after the Chi-Chi earthquake. At a trenching site in the southern segment as displayed with a cross in Fig. 1, Chen et al. (2004b) identified three pre-historic events (denoted by N-1, N-2, N-3). The 1999 event is denoted by N-0. Wang (2005) plotted the history of the four events and the cumulative slip (in meters) versus time (in years) for the four events. Results are displayed in Fig. 4 with the thick solid line segments representing the history of observed data. From the figure, the recurrence times are ~400 years for the 1999 event and ~850 years for the next one.

Stress Orientation Changes
From P-wave first motion polarities, Wu et al. (2010) determined the focal mechanisms of 4761 events during 1991 to 2007 in the Taiwan region. They found that the spatial distributions of stress axes were mainly controlled by tectonic structures; while the temporal variations were greatly influenced by the 1999 Chi-Chi earthquake. The orientation of the maximum horizontal compressive stress axes (SH) is generally in agreement with the direction of plate motion in a depth range 0 -30 km. They assumed that the trends of SH in the depth of 0 -10 km were strongly affected by the coseismic stress change of the 1999 Chi-Chi earthquake. In the northern half of the Chi-Chi rupture area, the trends of SH rotate 30° clockwise and the stress ratio increased by a factor of six after the mainshock. The orientations of SH still differ by 30° in 2007 comparing to that in the period before the Chi-Chi earthquake. The variation of SH is more diverse in the southern half of the rupture area, showing 20° counterclockwise rotation immediately after the mainshock followed by a clockwise rotation. The trend of SH returns to the pre-seismic direction of 110° in 2001. These notable changes of SH before and after the mainshock suggests that the magnitude of background stress in the rupture area is close to the coseismic stress drop. They also recognized a 10° counterclockwise rotation of SH in the entire ruptured area between 1991 and 1999 before the mainshock. However, to the south of the ruptured area, the trends of SH remained little changed before and after the mainshock. From the same data set, Hsu et al. (2010) also inferred that the spatial distribution of either the coefficient of friction or the pore pressure changed before and after the mainshock. From the two studies, the precursor time is ~9 years (listed in Table 1).

Temporal Variation in Seismicity Pattern
From the earthquakes with focal depths ≤ 20 km and M ≥ 3.4 reported by the Central Weather Bureau (CWB), Chen et al. (2005a) applied the pattern informatics algorithm (Rundle et al. 2003; and cited references herein) to recognize the seismicity pattern preceding the mainshock. Theirs results show anomalous changes of seismic activation and quiescence in a time period of ~6 years before the mainshock.
From the earthquakes with M ≥ 2 recorded by the CWB, Wu and Chen (2007) studied the temporal variation of monthly numbers of events from 1994 to 2006. From the results, they divided the whole time period into four intervals: I (from 1994 to 1999), II (from 1999 to 2000), III (from 2000 to 2002), and IV (from 2002 to 2006). They also calculated the standard normal deviate Z-value (Meyer 1975) for different intervals. Their results reveal that the ar-eas with relatively high seismicity in eastern Taiwan from 1994 to 1998 became abnormally quiet before the mainshock; while in the area with relatively low seismicity from 1994 to the occurrence of the mainshock in central Taiwan became unusually active after the mainshock. Wu and Chen (2007) assumed such a spatially changing seismicity pattern to be a phenomenon of 'seismic reversal' before and after the mainshock. Low seismicity in central Taiwan appeared at least 6 years before the mainshock. Thus, the precursor time for seismicity change is ~6 years (listed in Table 1).

b-Value Anomaly
Based on M ≥ 2.0 earthquakes with focal depths ≤ 40 km, Tsai et al. (2006) studied the temporal variations of b-values from 1 January 1994 to 31 August 1999 in three regions surrounding the source area (SA): region A including a small northern part of SA and a part to the north of SA, region B covering the main part of SA, and region C including a small southern part of SA and a part to the south of SA. They found that the highest b-value in region A appeared in early 1997 and the b-value in region B reached its peak in late 1997. Although there were relatively high b-values in region C during 1997, the b-value dramatically decreased to 0.65 after the M 6.2 Ruey-Li earthquake of 1 July 1998. They considered that the b-value anomaly in region C is likely contaminated by aftershocks of the Ruey-Li earthquake. Hence, they assumed that the b-value anomalies in regions A and B are a possible precursor of the mainshock. The precursor time for anomalous b-values is ~6 years (listed in Table 1).
From the catalog of CWB, Wu and Chiao (2006) selected 66069 M ≥ 2.0 earthquakes with focal depths ≤ 40 km from January 1994 to the mainshock. The data set that includes inland events and offshore ones located within 20 km of the coastlines of Taiwan is called the W&C date set hereafter. They calculated the b-values in the region before and after the mainshock. Results exhibit that the mainshock was preceded by a consistently decreasing trend of b-value as well as notable increases of seismic activity in regions surrounding the source area. The anomalous period started in January 1999 and lasted about 9 months, up to the occurrence of the mainshock. Since they only gave the time interval of decreasing b-value and did not provide the whole time period of anomalous b-values, the time length of 9 months is considered as the precursor time for b-value anomalies. Lee and Tsai (2004) measured the variations of P-wave travel-time residuals, δt P , before and after the mainshock by using seismic data from 1991 to 2002 reported by the CWB. Results show that the residuals increased at the stations immediately west of the Chelungpu fault about six years before the earthquake. The contour map of variations of residuals at station NSY near northern part of the Chelungpu fault reveals that most of P-wave ray paths passing through an anomalous zone east of the fault had positive residuals from 1994 to the mainshock. The anomalous zone is bounded by stations within 40 km to the west of and in the foot wall block of the Chelungpu fault. This implies that P-wave velocity, v p , began to decrease in such a zone in 1994, about six years before the mainshock. The precursor time for anomalous δt P or v p is ~6 years (listed in Table 1). Yu et al. (2001) observed the preseismic deformations of the earthquake in central Taiwan from annually repeated GPS data during 1992 to 1999. The total WNW-ESE shortening rate in the vicinity of the source area, ranging from the west coast to the western boundary of the Central Range, was up to 25 mm yr -1 . The pre-seismic crustal deformations were essentially a uniaxial compressional strain of 0.36 μstrain yr -1 along the direction of 114°. But, they did not provide the precursor time. Tsai et al. (2006) applied the ERS-2 radar images to detect the pre-seismic surface deformations in the surrounding area of the Chelungpu fault. They found that surface deformation began at least three years before the earthquake in an area immediately to the west of the northern segment of the Chelungpu fault where co-seismic surface deformations were also observed. The precursor time for surface deformations is ~3 years (listed in Table 1).

Seismic Quiescence
From the W&C data set, Wu and Chiao (2006) calculated the standard normal deviate Z-value (Meyer 1975) in the region before and after the mainshock. Results exhibit that the mainshock was preceded by a notable decrease of the regional seismicity rate. The anomalous period started in January 1999 and lasted about 9 months, up to the occurrence of the mainshock. The mean seismicity rate is 435 events per month with a standard deviation of 78 events. During the anomalous period, seismicity rates fell outside the range of one standard deviation with a mean value of 314 events per month. To clarify the problem of seismicity reduction, they also examined the monthly rate for events with M ≥ 4 only. But, they did not find significant reduction for such events in the period from 1 January to 19 September 1999. This indicates that the appearance of seismic qui-escence was attributed essentially to a remarkable decrease in the events with M < 4.
From the catalogue of the CWB, Wu and Chen (2007) observed that the monthly number of events remarkably decreased from January 1999 to the occurrence of the mainshock and then rapidly increased due to the appearance of aftershocks. From the spatial distributions of seismicity, the decrease in monthly numbers is particularly remarkable in the source area of the mainshock. Their observation is similar to that done by Wu and Chiao (2006). From the studies of the two groups, the precursor time for seismic quiescence is ~9 months (listed in Table 1).

Groundwater Level Changes
In 1999, there were 188 monitoring wells constructed on the Choshui River alluvial fan to the west of the Chelungpu fault by Water Resources Agency, Ministry of Economics, ROC (Chang et al. 2006). Chen et al. (2013Chen et al. ( , 2015 studied anomalous frequency characteristics of groundwater level before the earthquake from the corrected data after removing the factors of barometric pressure, earth tides, precipitation, and artificial pumping. They found that groundwater levels at 78% (= 42/54) of wells abnormally changed before the mainshock. For example, at the Huatang (HT) station 25 km to the west of the Chelungpu fault, groundwater level decreased about 250 days (about 17 January), reached the bottom (~-1.5 m) then immediately increased about 130 days (about 18 May), and returned to a local maximum 13 days (about 7 September) before the mainshock. At the HT station (shown with an open triangle in Fig. 1), the temporal variation in groundwater level from 1 August 1997 to 19 September 1999 is schematically displayed by thick line segments in Fig. 5 that is simplified from a figure in Chen et al. (2015). The dashed line denotes the regular decrease in groundwater level. For all stations in consideration, they observed that groundwater level changes ranged from 2 to 4 m. Their analyzed results show that amplitudes at the frequency band between 0.02/day and 0.04/day generally retained at the low stage and were apparently enhanced a few weeks before the mainshock. The precursor time for groundwater level changes is ~250 days (listed in Table 1).

Electromagnetic Precursors: Geomagnetic Annual Changing Rate
Since 1988, a geomagnetic network consisting of twenty-two repeated stations has been constructed in Taiwan by the Institute of Earth Sciences (IES), Academia Sinica (Yeh et al. 1981). Eight of them have also been equipped with a continuous recording system in which there is a proton precession magnetometer (i.e., Geometrics Model G-856) having a 0.1 nT sensitivity. The eight stations are shown with open circles in Fig. 1. Except for the Lunping (LP) station (shown by a solid square in Fig. 1) that is located about 130 km from the epicenter, others (shown by open squares in Fig. 1) are all located in the areas having high seismicity or remarkable crustal deformations. Other seven stations are: Liyutan (LY) and Tsengwen (TW) on western Taiwan; Neicheng (NC), Hualien (HL), Yuli (YL), Taitung (TT), and Hengchun (HC) on eastern and southern Taiwan. Chen et al. (2004a) assumed that the extra stresses acted on the source area before the mainshock could not produce significant effects on geomagnetic field at the LP station and thus they took this station as a reference one for others. They examined the time variations in the total geomagnetic field recorded at the eight stations from 1999 to 2001. They found that a zero isoporic zone (ZIZ), which is defined as the annual change rate of geomagnetic parameters ≤ ±5 nT yr -1 , appeared near the source area about 2 years before the mainshock. The precursor time for anomalous annual changing rate of geomagnetic field is ~2 years (listed in Table 1).

Groundwater Chemical Anomalies
After the Chi-Chi earthquake, Song et al. (2003) collected commercial bottled water (namely Chingjing water) that was pumped from wells at Puli, Nanton County in central Taiwan. Puli that is displayed with an open circle in Fig. 1 is about 18 km to the northeast of the epicenter. They analyzed the anions, e.g., Cl -1 , SO 4 -2 , and NO 3 -1 , of the bottled water in a time period of about 2.5 years from 1 December 1998 until after the mainshock. Results show that the concentrations of both sulfate (SO 4 -2 ) and nitrate (NO 3 -1 ) increased steadily from 1 April 1999, as compared to the average constant level measured from 1 December 1998 to 31 March 1999. The concentrations reached an excess of 129.9 and 94.7% in April 1999, remarkably dropped in July 1999, and continued to decrease until the mainshock. The concentrations decreased to about 10% for sulfate and 20% for nitrate that were below the respective constant values after the mainshock. However, they did not report the variation in Cl -1 concentration. The precursor time for groundwater chemistry anomalies is ~7 months (listed in Table 1). Wang et al. (2005) analyzed some isotopic and hydrological changes related to the earthquake. They reported the changes in the isotopic composition of the Choshui River alluvial fan near the Chelungpu fault and the groundwater level in the fan during and following the mainshock. They reported three aspects of hydrological changes. First, following the earthquake the lower aquifers beneath the fan showed a significant shift in isotopic composition towards that of the surface water in the river, suggesting enhanced exchanges of water between the river and the ground water. Second, in some wells water levels and isotopic compositions in different aquifers converged to the same respective values during the earthquake, suggesting coseismic exchanges of water between the different aquifers, which implies enhanced permeability due perhaps to the fracturing and breaching of aquitards between the aquifers. Third, the pattern of the coseismic water-level response is distinctly different from that of the shift in the isotopic composition, suggesting that they were produced by different mechanisms. Although the authors found some isotopic and hydrological changes, they only reported the average values in a long time period and did not provide the most possible times of changes. Hence, it is difficult to estimate the precursor time of isotopic and hydrological changes from their study.  Yen et al. (2004) analyzed the data recorded at eight geomagnetic stations of the IES (see Fig. 1). Although a magnetic storm occurred two days after the mainshock and was recorded at the LP station, not any disturbance in the total geomagnetic intensity (TGI) was found at the LP station before and during the mainshock. Their results reveal that there were significant fluctuations, with the largest amplitude up to 200 nTs, in the differences of TGI between the LY station, whose epicentral distance is shorter than 50 km, and the LP station from mid-August (about 1.1 months before the mainshock) to November 1999 (about 2 months after the mainshock). The precursor time for geomagnetic anomalies is ~1.1 months (listed in Table 1).  observed seismo-geomagnetic anomalies before the mainshock at the LY station during 1988 to 2001. They computed the diurnal range ratio (DRR) between the LP station and the LY station from the data. They also calculated the ratio of monitored number for each DRR to the total monitored number in five different time periods before and after the mainshock and the average ratio of monitored number for each DRR in the whole thirteen years. The average ratio is taken to be the reference. Finally, they plotted the distributions of the ratio in five different time periods and the reference. Results show that the distribution in the time period just one month before and one month during and after the mainshock significantly departs from the reference distribution. To interpret the departure, they assumed that changes of underground conductivities and currents around and below the epicenter significantly affected the nearby preseismic ground geomagnetic field. Of course, such a distribution is highly dependent upon the focal mechanism. In comparison with the results obtained by Yen et al. (2004), the precursor time for seismo-geomagnetic anomalies is ~1 month. Gokhberg et al. (1982) first took the EM emissions in the low-frequency (LF) band (30 -300 kHz) as a precursor of earthquakes. Since then, EM signatures in different frequency bands, i.e., from extremely low frequency (ELF) band for 3 -30 Hz to very high frequency (VHF) band for 30 -300 MHz, have been considered as a promising candidate of a short-term precursor because of convincing reports about the appearance of related signatures before several large earthquakes (Ohta et al. 2001; and cited references therein). Among them, the signatures in two lower frequency bands, i.e., ELF and VLF (3 -30 kHz), are particularly promising. For the Chi-Chi earthquake, Akinaga et al. (2001) measured the polarization that is the ratio, Z/G, of vertical magnetic field, Z, to the horizontal one, G, from ULF (300 -3000 Hz) emissions recorded at the LP station. They found that the Z/G significantly increased about two months before the mainshock, and thus they considered this change as a precursor of the earthquake. The precursor time for anomalous ULF signatures is ~2 months (listed in Table 1).

Slow-Slip Events
Lin (2012) measured the surface crustal deformations from the integration of broadband velocity seismograms of the CWB. From the results, he found significant deviations of the vertical displacement from a normal Earth tidal pattern during 15 to 19 September before the mainshock. Meanwhile, he also observed a series of slow-slip events occurring on the nearly horizontal plane (i.e., the decollement) of the Chelungpu fault at depths between 10 and 12 km. Although the hypocenters of these slow-slip events are not well constrained owing to limited observations only at two seismic stations, he assumed that these events played the primary role on producing anomalous surface crustal deformations and could be a precursor. The precursor time for slow slip is ~5 days (listed in Table 1).

Infrasound
In Institute of Acoustics, Chinese Academy of Science at Beijing, PRC, there is an infrasound recording system that is located at (39.87°N, 116.48°E) and consists of a CC-1T type capacitive sensor with power supply and data logger. The sensor has a sensitivity of 0.01 Pa in the frequency range 0.5 -200 Hz and a dynamic range of 80 dB (Xie 1991). In the system, the analog and digital records have been continuously operated since 1982 and 2002, respectively. Xia et al. (2011) observed anomalous infrasound signals before 92 M ≥ 7.0 worldwide earthquakes, including the Chi-Chi earthquake, during 2002 to 2008. For the Chi-Chi earthquake, they observed a peak amplitude of 1100 mV, which is equivalent to a sound pressure of 8.8 Pa, from the recorded signals at 16:00 -16:40 pm (Beijing local time) on 18 September (about 3 days) before the mainshock. The precursor time for infrasound signal is ~3 days (listed in Table 1). Liu et al. (2001Liu et al. ( , 2004a analyzed the seismo-ionospheric signatures prior to the mainshock by examining the TEC recorded by the receivers of a network of the global positioning system (GPS) in Taiwan. They also analyzed the temporal variations in TEC recorded by a sweep frequency pulsed radar device, i.e., ionosonde, located at Chungli (CL) (25.0°N, 121.l°E) whose location is almost the same as the LP station. Chungli is shown with an open circle in Fig. 1. There is a similar tendency between the two types of observations. They combined the data of 13 GPS receivers together to examine the temporal and spatial variations of TEC prior to the mainshock. Results show that the equatorial anomaly crest moves equator-ward and the significantly decreased before the mainshock. Based on Liu et al. (2001Liu et al. ( , 2004a, a temporal variation in TEC is schematically displayed in Fig. 6a in which the solid line and dotted lines represent the observations and references (previous 15-day median), respectively. This figure is simplified because the fluctuations are not included. A comparison between the disturbed and reference data suggests that TEC decreased significantly in the afternoons about 4, 3, and 1 days before the mainshock.

TEC and f o F 2 Anomalies
The virtual height of the ionosphere may be measured from remote sensing using radio waves recorded by an ionosonde. The virtual height is equivalent to the product of onehalf the time-of-flight of the transmitted radio wave and the speed of light. The plot of frequency versus virtual height is called an ionogram that displays the frequency-dependent variation of the virtual height of reflection. Usually, there are two traces, i.e., O-mode and X-mode, on an ionogram. Based on the magneto-ionic theory (Budden 1985), the plasma frequency is equal to the vertically reflected O-mode frequency. The highest (or critical) frequency, f o F 2 , on an O-mode trace may be considered as the penetration plasma frequency or the largest density of the ionosphere. Chuo et al. (2002) analyzed the f o F 2 observed by the CL ionosonde. By using the 15-day running mean and the associated standard deviation computed before the mainshock, they estimated the upper and lower bounds and then detected the seismo-ionospheric perturbations. Their results are similar to Fig. 6a. Results show that the perturbation appeared few days before the earthquake. From the two studies, the precursor time for TEC and f o F 2 anomalies is 3 -4 days (listed in Table 1).

Lightning: An Atmospheric Anomaly
From the records of cloud-to-ground lightning occurred 15 days before and after the mainshock, Tsai et al. (2006) found that the frequency of lightning significantly increased on 17 September 1999, which is 4 days before the mainshock, and on this day the lightning occurred mainly near the southern end of the Chelungpu fault and thus near the mainshock epicenter. The precursor time for anomalous lightning is ~4 days (listed in Table 1). Ohta et al. (2001) analyzed the ELF/ULF emissions recorded at Nakatsugawa observatory, that is located at (35.4°N, 137.5°E) in Gifu Prefecture, Japan, from 1 January 1999 to 21 September. The ELF noise level increased by more than 5 dB from the normal level in 1.5 hours during 21:30 -23:00 pm Japanese Standard Time (or 20:30 -22:00 pm Taiwan Local Time) on 20 September and the upper limit extended up to 50 Hz. Through a careful comparison of their observation with the nearby lightning as detected by VLF, they confirmed that the abnormal increase in ELF noise level was not due to the lightning. The measured phase difference of the ELF emissions (B X , B Y ) indicates that these emissions were linearly polarized and thus had propagated in the sub-ionospheric waveguide over a long distance. From the identification of the goniometric direction from the polarization, they indicated that the main direction of the ELF/ULF emissions was pointing toward Taiwan. This made them to assume that the ELF/ULF emissions were produced from the preseismic slip of the Chelungpu fault before the mainshock. The precursor time for anomalous ELF/ULF signals is about 4 hours (listed in Table 1).

Sky and Earthquake Lights
From interviews with local people, Chen et al. (2000) collected the data of preseismic sky light (with different colors) and coseismic seismic (green) light. The 'seismic light' is usually called the 'earthquake light' in the seismological community. The 'sky light' might also be the earthquake light because it was associated with the earthquake. Preseismic sky lights were seen by local people several times from north to south along the Chelungpu fault. The co-seismic earthquake light was seen by local people only once at a site in the northern segment of the fault. This difference will be discussed below. The precursor time for sky light is few hours (listed in Table 1). Song et al. (2006) collected water samples from both hot and artesian springs in Kuantzeling (displayed with an open circle in Fig. 1), Chiayi in west-central Taiwan from 15 July 1999 to the end of August 2001 and measured cation and anion concentrations. Results reveal that the concentrations of chloride and sulfate ion abruptly increased on 19 September about two days prior to the mainshock and lasted a few days afterward. These anomalies are characterized by remarkable increases in Clconcentrations above the means at both hot and artesian springs. However, the SO 4 2concentrations did not change on 19 September. The precursor time for anomalous Clconcentrations is ~2 days (listed in Table 1).

Biological Precursors: Anomalous Animal Activities
From interviews with local people, Chen et al. (2000) collected the data of anomalous activities for 12 kinds of animals at 28 locations before the mainshock. Except for the 28 th location at Jiou-Fen-Erl-Shan with an epicentral distance > 10 km, other 27 locations are very close the Chelungpu fault. Jiou-Fen-Erl-Shan is almost in between the Chelungpu fault and Puli. The aberrant behavior of ants occurred as early as 8 -10 weeks at a location and 3 days at four other places before the mainshock. The aberrant behavior of cicadas occurred 4 -6 weeks at a locality before the mainshock. The aberrant behavior of earthwarms, diplopods, and fishes occurred about 1 -2 weeks at some locations before the mainshock. The aberrant behavior of birds occurred ≤ 7 days at three locations before the event. The roachs abnormally appeared 3 days at a location before the mainshock. The cats abruptly disappeared and turtles abruptly appeared at the same local area ~1 day before the mainshock. The aberrant behavior of dogs occurred at several locations < 1 day before the mainshock. The snakes abruptly appeared at a location ~2 hours before the mainshock. The precursor times for all animals in study are listed in Table 2 where there are three time periods: week, days, and hours.

Some Pre-Seismic Natural Phenomena
From interviews with local people, Chen et al. (2000) collected the data of some anomalous natural phenomena, including wind, sound, smell, and initial motions at 28 locations before and during the mainshock. They divided the data into two types: (1) pre-seismic phenomena: wind; and (2) co-seismic phenomena: sound, smell, and initial motions. The changes of wind directions occurred on 20 September, i.e., the day immediately before the mainshock. Since several factors, including the atmospheric pressures, temperature, raining etc., may control the directions of wind, I wonder if the changes of wind directions may be considered as a precursor or not.

Very-Long-Term Prediction: Earthquake Recurrence
As shown in Fig. 2, earthquake recurrence shows (regular or irregular) repetitive occurrences of earthquakes on a fault. As shown in Fig. 3, the strain energy almost totally released during the last event, then accumulated again under increasing stress loading starting at time t 0 in Stage 1, and finally release again in Stage 4 starting at time t r . The recurrence time is from several ten to hundred years, even several thousand years. In a tectonically active region, a fault is exerted by relative movement of two plates with a speed of V p as displayed in Fig. 2. Assume that V p be applied at a distance w from each side of the fault, the strain is ε = V p t/2w. Thus, 2w is the width in which the shear strain develops progressively across the fault prior to an earthquake. The static stress, σ s , is σ s = σ d + μV p t/2w, where σ d and μ are the dynamic friction stress and the crustal rigidity, respectively. When σ = σ s at T R = 2wΔσ/μV p , where Δσ = σ s -σ d is the stress drop, the fault slips. Hence, T R is the inter-event time or the recurrence period (or time) of an earthquake sequence. Shimazaki and Nakata (1980) proposed three simple earthquake recurrence models (see Fig. 7), each with a constantly increasing tectonic stress. The three models are: (1) the perfectly periodic model; (2) the time-predictable model; and (3) the slip-predictable model. The detailed description and debates about the three models can be found in Wang (2019). Numerous factors, including wear process, friction, plasticity, strain rate etc., can influence the pattern of earthquake occurrence. Wear process (Wang 2018b) and time-strengthening static friction (Wang 2020) on a fault are two major factors in influencing the pattern of earthquake recurrence.
In order to meet the slip-predictable model represented by a thin dotted line in Fig. 4, it is necessary to move the history of slip with a time shift to the left as displayed by the thick dashed line segments in the figure. It sounds unreasonable because the occurrence times of historical events are incorrect. The lower thin solid line and the upper dotted line display the long-term slip rates of V s = 4.97 mm yr -1 for the time-predictable and slip-predictable models, respectively. The minus sign denotes "BP." The values of T R of N-0 and N-3 are two controlling points of the plot. From the time-predictable model, the optimum values of T R are about 400 years BP for N-1 and 980 years BP for N-2. The estimated T R of N-1 that is within the range of 430 -150 years BP obtained by Chen et al. (2004b) or within that of 500 -300 years BP inferred by Ota et al. (2001) is acceptable. Chen et al. (2004b) expressed that N-2 might occur in 700 -800 years BP, with a higher uncertainty. The value of 980 years BP seems to be a good one. Hence, Wang (2005) assumed that the time-predictable model is acceptable. The next earthquake will happen 850 years in the future. The study area has been steadily deforming since 0.7 Ma ago. About 30% strains caused by regional tectonic loading were released during failures of the Chelungpu fault. From the data obtained from other trenching sites, Chen et al. (2005b) also assumed the time-predictable earthquake recurrence along the fault.
Based on a one-body model, with the surface displacements estimated at several near-fault sites, Wang (2003) proposed that the recurrence period of the northern segment of the fault is about 3 times longer than that of the southern one. This suggests that the next Chi-Chi earthquake will rupture only along the southern segment of the Chelungpu fault. However, he could not evaluate the definite recurrence period from the model. Since the recurrence time is very long, the earthquake recurrence is more appropriate for earthquake forecasting than for earthquake prediction.

Long-Term Prediction
As shown in Fig. 3, long-term precursors started at t 0 and may last for a long time including Stages 2 and 3. During a long time period, the fault zone is under elastic loading in Stage 1 and plastic strain hardening in Stages 2 and 3 due to dilatancy, which is an increase in volume prior to failure (e.g., Nur 1972;White 1976), in the last two stages because of crack coalescence and fluid transport in the fault zone. The related precursors are discussed below. change from higher seismicity to lower seismicity in Taiwan over a time period of ~6 years before the mainshock. They also saw two significant phenomena as follows: First, to the east of the Chelungpu fault seismicity changed from a relatively high state during1994 to 1998 to a relatively low one in the first nine months of 1999 before the mainshock. Secondly, to the west of the Chelungpu fault seismicity changed from a relatively low state from 1994 to the occurrence of the mainshock to an unusually active one after the mainshock. They called such a change 'seismic reversal' with the changing point at the occurrence time of the mainshock. 'Seismic reversal' was first qualitatively defined by Russian seismologists (Keilis-Borok et al. 1994;Shebalin et al. 1996) as follows: zones of relatively high seismicity become unusually quiescent or zones of relatively low seismicity become unusually active. This commonly takes place a few months prior to an impending large earthquake within a distance of about 100 km from its future epicenter. Several authors (e.g., Scholz 1988Scholz , 1990Main and Meredith 1991;Knopoff 1996b) proposed different physical models to interpret the phenomenon. To the east of the Chelunpu fault, the change was mainly due to a decrease in seismicity just around the source area; while to the west of the Chelunpu fault, the change was mainly due to an increase in aftershocks. Hence, I wonder if the change of seismicity pattern as observed by Chen et al. (2005a) and Wu and Chen (2007) is indeed seismic reversal or not, because an increase in seismicity after the mainshock is due to the occurrences of aftershocks.

Variation in Seismicity Pattern: Seismic Reversal and Seismic Quiescence
As mentioned above, Wu and Chiao (2006) and Wu and Chen (2007) found a seismic quiescence started in January 1999 and lasted about 9 months until the mainshock. Professor Mogi of Tokyo University, Japan first studied seismicity patterns, including seismic quiescence, prior to earthquakes with M ≥ 6 (Mogi 1981; and cited references therein). Since then, many studies on seismic quiescence before small to large earthquakes have been reported. In fact, this kind of study has been lasted until now (e.g., Kanamori 1981;Habermann 1988;Chen et al. 1990;Shebalin et al. 1996;Shebalin and Keilis-Borok 1999;Zöller et al. 2002;Sukrungsri and Pailoplee 2017;Peresan 2018). Essentially, there are three types of seismic quiescence (Scholz 1988;Main and Meredith 1991): intermediate-term, short-term, and post-seismic seismic quiescence. Habermann (1988) reviewed the studies on seismic quiescence done before 1988. Clearly, the seismic quiescence identified by Wu and Chiao (2006) and Wu and Chen (2007) is the intermediate-term one. Previous studies reveal that seismic quiescence is due to a decrease in the events with magnitudes smaller than the magnitude of the impending earthquake. But, seismic quiescence claimed by Wu and Chiao (2006) was caused mainly by a remarkable decrease of micro-events with M < 4. The occurrences of micro-events with M < 4 are usually not felt by most human beings, except for some of those living in the source area, and behave almost quiescent before or after the mainshock. Although this phenomenon is not typical, it is very interesting. The reason to cause a remarkable decrease in micro-events with M < 4 before the mainshock is still open, and thus a mechanism should be developed to interpret the phenomenon.

b-Value Anomaly
The temporal variation in b-value is considered as one of significant precursors for volcano activities and earthquake occurrences. Gorshkov first observed a decrease in bvalue before the Russian Bezymianny volcanic eruption in 1956 (see Aki 1985). Suyehiro (1966) first found a change of the b-value before and after an earthquake and thus considered the b-value anomaly to be a precursor. This phenomenon was also observed by numerous researchers (see Chen et al. 1990;Wang et al. 2015Wang et al. , 2016; and cited references therein). The precursor time of b-value anomaly inferred by Tsai et al. (2006) is ~6 years. As displayed in Fig. 3 For the Chi-Chi earthquake with M = 7.6, the value of T is ~4 years. Chan et al. (2012) also measured the average value of T (≈3 years) for 23 Taiwan's earthquakes with M L ≥ 6. Although the Chi-Chi earthquake was included in their data set, they did not give the value of T. Their average T is slightly shorter than that calculated from Eq. (1) and much shorter than that given by Tsai et al. (2006). Since the data points for larger earthquakes in Wang et al. (2016) are somewhat scattered, the estimated value of 6 years by Tsai et al. (2006) is still in the range of standard errors and thus acceptable.

Temporal Changes in P-Wave Travel-time Residuals
Based on seismic-wave travel time data, Russian seismologist Semenov (1969) first claimed that P-and S-wave velocities (denoted by v p and v s , respectively, hereafter) and related quantities decreased by 10 -20% before an earthquake and then recovered. He also stressed that this recovery followed by an impending earthquake. Wyss (1975) reported this precursor from Russian to English. Later, numerous researchers (see Geller 1996; and cited references therein) studied anomalous P-and S-wave travel-time residuals (denoted by δt p and δt s , respectively, hereafter). For example, Robinson et al. (1974) observed a sharp increase of v p about two months before an M 5 earthquake and they suggested that the increase was due to dilatancy. However, Lindh et al. (1978) claimed that this was an artefact. Whitcomb (1976) reported a velocity anomaly beneath the Transverse Ranges in southern California. This anomaly was considered by Hammond (1976) and Shapley (1976) as a precursor of an impending event with M = 5.5 -6.5 in the area. But, no such event happened. Some authors (e.g., Nur 1972;Scholz et al. 1973;Whitcomb et al. 1973;Anderson and Whitcomb 1975;Griggs et al. 1975) proposed that dilatancy could explain the decrease in v p and v s .
From controlled sources (quarry explosions), several authors obtained negative results. Allen and Helmberger (1973) found no significant variations in v p before an earthquake. Johnson (1973, 1974) found temporal variations of at most 1%, which could be explained without invoking in situ changes. Kanamori and Fuis (1976) found temporal variations of the order of at most 1%. Some studies (e.g., Leary and Malin 1982;Haase et al. 1995) claimed only smaller bounds on temporal variations. Aggarwal et al. (1975) reported late arrivals (by 0.13 sec, corresponding to 0.04 mm recorded on the smoked paper) from quarry explosions concurrent with a low in t s /t p in earthquake data. This may be at the noise threshold. Noted that since ~1980, only very few reports of large velocity changes have been published. Lee and Tsai (2004) reported an increase of δt p or a decrease of v p in an area to the west of the Chelungpu fault about six years before the earthquake. In previous studies, the precursor times for changes in v p and v s are usually only several months. But, the precursor time inferred by Lee and Tsai (2004) is ~6 years and similar to that of b-values estimated by Tsai et al. (2006). Wang (2016) proposed a mechanism to interpret the temporal variation in b-value before an earthquake. A detailed explanation can see Wang (2016) and thus only a simple description is given below. Laboratory experiments (Cadoret et al. 1995) exhibit that the v p of rocks is strongly controlled by the water saturation in the rocks. In Fig. 6, the solid line schematically shows the general pattern of variation of v p in terms of water saturation as well as time. Based on the fault-valve model proposed by Sibson (1992), the water inside a fault zone almost completely released and the zone is sealed after an earthquake. At t = t 0 (see Fig. 6), the water is re-ejected into the fault zone from the upper mantle or/and from the crust. This process works on both the main fault (for the mainshock) and many subfaults (for forerunners and foreshocks). Since the permeability of rocks inside the source area is commonly lower than 10 -12 m 2 (Wang et al. 2009), the process of water ejection into the fault zone is very slow. The value of t 0 from last event could be tens or hundreds or thousands of years and depends on the size of the magnitude of the earthquake occurring on the fault. Water re-ejection increases the water saturation in the fault zone. When t > t 0 , the amount of ejected water is high enough to substantially arise the pore pressures inside the voids of fault rocks. This makes the voids continuously open and speed up water ejection, thus remarkably rising degree of water saturation in the fault zone. This results in a decrease in v p and an increase in P-wave travel time as well as δt p , in comparison with normal v p of dry rocks. At t = t m , water saturation is increased to a certain critical level and thus the voids cannot be opened any more. This leads to the minimum v p as well as the maximum δt p . When t > t m , v p increases and thus δt p decreases. At t = t v , v p returns to the normal value and thus δt p = 0. The time interval from t 0 to t v may be several years, depending on the size of impending earthquake. When t > t v , the degree of water saturation within the fault zone becomes very high and v p may be higher than that of dry rocks. This leads to δt p < 0. This process might occur in step 3c when fluid diffusion happens. Theoretically, Wang (1995) found that the b-value decreases with increasing elastic parameter of a fault-zone materials. This made Wang (2012) obtain a relationship between the elastic parameter and v p , From the two correlations, Wang (2016) related the temporal variation of b-value (shown by a dotted line) to that of v p as displayed in Fig. 8. Like v p , anomalous b-value started from the early time of Stage 2. The mainshock occurs at t = t r , and thus T = t rt 0 is the precursor time of anomalous b-value, v p , and δt p .
The observed temporal variation in b-value (Tsai et al. 2006) and that in δt p (Lee and Tsai 2004) seem able to be comparable, respectively, with the modeled ones in b-value and δt p , even though there are fluctuations in both observed curves. From Tsai et al. (2006) the occurrence of the peak b-value, especially for region B, appeared between 1997 and 1998; while from Lee and Tsai (2004) there were several peak values of residuals because of high fluctuations of data. The study areas from which the two groups of authors took the data are different. The study area for b-values is larger than and covers that for the residues. This might cause the difference in the patterns of temporal variations for the two parameters. Lee and Tsai (2004) also assumed that the anomalous zone resulting in the decreases in v p is bounded by stations within 40 km to the east of the Chelungpu fault and in the foot wall block of the fault. From wide-band magnetotelluric (MT) surveys, Chen et al. (2007) inferred a low-resistivity (LR) zone, beneath the hypocenter and to the west of the Chelunpu fault. They assumed that this zone was the source of deep-crustal fluids which might be one of the factors in producing the rupture processes of the Chi-Chi earthquake. The depth range of this LR zone is almost similar to the anomalous zone claimed by Lee and Tsai (2004). Sibson (1977Sibson ( , 1992 mentioned that water re-ejection may happen in the whole fault zone rather than in a smaller one. The anomalous zone claimed by Lee and Tsai (2004) could be the main one to affect δt p , yet not the only one.

Intermediate-Term Prediction
As shown in Fig. 3, intermediate-term precursors appeared at t y and may last for a long time from Stages 2 to 3. The related precursors are discussed below.

Crustal and Surface Deformations
The preseismic deformations in central Taiwan observed by Yu et al. (2001) do not show any remarkable change of deformations before the mainshock. While, the preseismic surface deformations in central Taiwan observed by Tsai et al. (2006) exhibit that surface deformation began at least three years before the earthquake in the area immediately to the west of the northern segment of the Chelungpu fault. From the south to the north along the fault, the measured co-seismic surface displacements range from 1.0 to 11.1 m for the horizontal component and from 2 to 7.5 m for the vertical component (cf. Wang 2006). It seems to suggest the existence of a correlation between preseismic slip and co-seismic one. From Fig. 3, the surface deformations might start from the beginning or the early time of Stage 2 when the stresses in the source area were high enough to yield such deformations. In other word, Stage 2 started about 3 years before the mainshock. Chen et al. (2013Chen et al. ( , 2015 observed unusual fall and rise of groundwater levels at 78% of wells about 250 days (> 8 weeks) before the mainshock. From the temporal variation in groundwater level (a simplified example as displayed in Fig. 5), they defined two phases in the time interval having the fall-and-rise process of groundwater level: Phase 1 for the fall of level from the beginning of a decrease to the bottom and Phase 2 for the rise of level from the bottom back to the normal value.

Groundwater Level Changes
The large precursory water level changes have been observed in confined aquifers (e.g., Roeloffs et al. 1997). For such aquifers, Rice and Cleary (1976) related the change in reservoir fluid pressure, Δp, to the incremental change in volumetric strain, Δε, that is positive for tension (or dilatation) and negative for compression in the following equation: where μ is the shear modulus, C s is Skempton's coefficient, and ν u is the undrained Poisson's ratio. The change in water level Δh is related to Δp by where t is the fluid density and g is the gravitational acceleration. Hence, a fall in water level, i.e., Δh < 0, related to tensional strain (Δε > 0) because of tensional pressure (Δp < 0); while a rise in water level, i.e., Δh > 0, is corresponding to compressive strain (Δε < 0) due to a compressive pressure (Δp > 0). For typical values of μ = 3 GPa, C s = 0.8, and ν u = 0.3, the water level change is 52 cm per 10 -6 strain (Roeloffs 1988). On the other hand, for unconfined aquifers, the water level change is given by where H s is the saturation thickness of the aquifer and } is the porosity. For a 100 m saturated aquifer with } = 0.02, the expected change in water level is 0.5 cm per 10 -6 strain (Roeloffs 1988). It is much smaller than that for a confined aquifer. Chen et al. (2013Chen et al. ( , 2015 mentioned that unusual fall and rise of groundwater levels occurred at open wells. Hence, the water level changes happened in unconfined aquifers as described by Eq. (4) and was corresponding to tensional strain in Phase 1 and to compressive strain in Phase 2. However, the Chi-Chi earthquake was caused by thrust faulting with a strong strike-slip component because central Taiwan is under compressive environment due to regional tectonics. Chen et al. (2015) measured the daily strain rate, ε', from the GPS data. They defined ε' < 0 for the compressive stress and ε' > 0 for the tensional stress. They also inferred the time-varying spatial distributions of stresses in the region surrounding the Chelunpu fault from the GPS data. From a comparison between the spatial distribution of groundwater levels and that of ε', they assumed that in Phase 1 the tensional stress (with ε' > 0), which was yielded in the environment of a large regional compressive stress, resulted in a fall of groundwater level; while in Phase 2 the compressive stress (with ε' < 0) led to a rise in groundwater level. From Fig. 3, the groundwater level changes might start from the later time of Stage 2 when the stress in the source area is approaching its peak. Song et al. (2003) observed that groundwater chemistry anomalies appeared about 7 months before and retained after the mainshock. They assumed that the temporal variations in SO 4 2and NO 3were not due to storage conditions, including dissolution, absorption, or evaporation of water in the bottles, and due to other reasons. Several natural factors in influencing changes in the chemical compositions of groundwater are: (1) including different compositions of groundwater recharge, the petrologic and mineralogical compositions of subsurface rocks, and water-rock interactions (Domenico and Schwartz 1990;Langmuir 1997); (2) fluid-source switching in response to fault sealing and unsealing, with the newly tapped aquifer containing chemically and isotopically distinct water and the changing stress state associated with the mainshock (Claesson et al. 2004); (3) the mixing of different water compositions (King et al. 1981;Thomas 1988); and (4) stress-induced electro-chemistry (Paudel et al. 2018). Of course, artificial pollutants may also be a factor. Domenico and Schwartz (1990) claimed that sulfate ion is unstable under groundwater conditions and may be changed by several reactions: (1) sulfide mineral oxidation, the precipitation-dissolution of gypsum in an unsaturated zone; and (2) the dissolution of anhydrite or gypsum or by redox reactions in a saturated zone. Although the nitrate ion is more stable than sulfate ion, but it can also be changed by redox reactions (Domenico and Schwartz 1990). Hence, Song et al. (2003) suggested that the temporal variations in concentrations of SO 4 2and NO 3are mainly caused by the mixing of water bodies with different compositions came from distinct sources, yet not due to artificial pollutants.

Ground-Water Chemistry Anomalies
These preseismic chemical changes may be attributed to the stress/strain-induced pressure changes in the subsurface water system, then followed by limited precursory geochemical discharges generated by limited changes in the levels of the subsurface reservoirs, finally led to the mixing of deeper aquifer. The groundwater was diluted by the injection of surface water into the subsurface system. This confirms the possibility by using groundwater chemistry as a precursor of an impending earthquake. From Table 1, we can see that the precursor time of groundwater chemistry anomalies is comparable with that of groundwater level changes. The Puli area is about 18 km away from the Chelungpu fault and thus almost in the deformed zone of the earthquake. The time-varying spatial distributions of daily strain rates inferred by Chen et al. (2015) shows that the strain rate in the Puli area became negative about 200 days before the mainshock. This means that the loading stresses on the underground rocks were compressive. The compressive stresses may decrease the porosity of affected strata, thus pushing the fluids in cracks flow into the aquifer. Upwelling of groundwater from the deeper rocks may carry the chemical components that originally existed at depths up to the shallow part, thus leading to the mixing of deeper aquifer and causing groundwater chemistry anomalies. About 130 days before the mainshock, the strain rate in the area became positive and thus tensional stresses loaded on the underground rocks. This increased the porosity of affected strata, thus pushing the fluids in aquifer back into the aquifer. This reduced the mixing phenomenon and geochemical anomalies decreased as observed by Song et al. (2003). From Fig. 3, the groundwater chemistry anomalies might start from the later time of Stage 2 when the regional stress was approaching its peak.

Short-Term Prediction
As shown in Fig. 3, short-term precursors occur mainly due to microcrack linkage and pore fluid diffusion almost in the steps 3a and 3b of Stage 3. The short-term precursors appeared in the later time of step 3a and the whole step 3b. The time window is from six months to eight days before the earthquake. Related precursors are discussed below.

Electromagnetic Precursors: Geomagnetic Anomalies and ULF Emissions
Preseismic EM emissions including ULF signatures and anomalous geomagnetic field have been considered to be caused by either direct radiation from a fault zone (e.g., Ogawa et al. 1985;Molchanov and Hayakawa 1995;Molchanov et al. 2004) or a change in geoelectric conductivity inside and near a fault zone. The latter leads to a change in ULF waves generated by magnetospheric sources (Merzer and Klemperer 1997). Several models were proposed to interpret the generation of EM emissions (see Molchanov and Hayakawa 1995). The details concerning observations and the models of EM emissions can be found in Shrivastava (2014). Since ULF emissions observed before and after an earthquake behave as transient variations, some authors (e.g., Hayakawa 1994, 1995;Molchanov et al. 2004) proposed that microfracturing electrification to be the most possible mechanism to generate such emissions. Using a space-time simple model of microfracturing, they compared theoretical results and observations especially for the intensity, spectrum, and temporal development of ULF magnetic field variations. From a positive result of comparison, they confirm the mechanism. Akinaga et al. (2001) observed the ULF signatures about 2 months before the Chi-Chi earthquake. Meanwhile, Yen et al. (2004) and  observed the anomalous geomagnetic field about 1 month before the mainshock. Tsai et al. (2006) emphasized the consistence between the two EM phenomena. Nevertheless, two problems arise. The first problem is: Why did the ULF appear one month earlier than the anomalous geomagnetic field, because the EM emissions are usually associated with the geomagnetic and electric fields?
The answer is open. The second problem is: From the temporal variation in daily amounts of precipitation from 1 December 1998 to 30 April 2001 in the Puli area of central Taiwan as displayed in Fig. 2c of Song et al. (2003), high precipitation appeared from 1 May to 30 September of 1999. This time period is almost the regular raining season in Taiwan. Although the Puli area is about 10 km to the east of the Chelungpu fault, the precipitation data recorded there may be considered as a representative of the raining season in central Taiwan. Raining is often accompanied with strong cloud-to-ground lightning. An example can be seen from Fig. 9 of Tsai et al. (2006). That figure exhibits the appearance of island-wide cloud-to-ground lightning on 17 September 1999 just three days before the mainshock. Lightning can produce EM anomalies. Hence, it is necessary to examine the possible effect of cloud-to-ground lightning on the disturbance of geomagnetic field and EM emissions. Only after removing such an effect, it is possible to directly correlate the two kinds of EM precursors. Hayakawa (2013) mentioned that in general ULF emissions first appeared with an intensity enhancement about 1 -2 weeks before an earthquake and then lasted for about one week (at least a few days) before an earthquake. The emissions increased again a few days before the earthquake and then were followed by an abrupt increase only a few hours before the earthquake. For the Chi-Chi earthquake, the ULF signatures appeared about 2 months before the event, lasted for more than one month, and finally disappeared a few days before the event (Akinaga et al. 2001). A time interval of 2 months (~8 weeks) is much longer than that of 1 -2 weeks. In addition, the UFL emissions were not re-active after quiescence. From Fig. 3, ULF emissions started from the earlier time of step 3a and continue to step 3c. Since the pre-seismic slip is small in steps 3a and 3b, the intensity of ULF emissions should be weak in this time period. In case that the emissions last to step 3c, the intensity must become stronger than before. Indeed, co-seismic ULF emissions may exist because co-seismic slip is much longer than pre-seismic slip. From the observations by Akinaga et al. (2001), the intensity of ULF emissions decreased and somewhat disappeared a few days before the earthquake. The reason to cause such a temporal variation is unclear.

Imminent Prediction
As shown in Fig. 3, imminent precursors occur due to pore fluid diffusion and quasi-static slip almost in the very later time of step 3b and the whole step 3c of Stage 3. The related precursors that appeared almost in the seven days before the mainshock are discussed below.

Slow-Slip Events
Lin (2012) observed the seven slow-slip events on the fault plane of the Chelungpu fault beneath central Taiwan about 5 days prior to the mainshock. He considered those events as a precursor of the Chi-Chi earthquake. Nevertheless, he also emphasized the importance of a basic condition that the fault plane should be the decollement from a geological standpoint. The evidence to confirm the fault plane to be the decollement beneath central Taiwan can be found from numerous studies (e.g., Suppe 1984;Kao and Chen 2000;Chen et al. 2002;Kim et al. 2010).
In his study, there is an interesting phenomenon that the inferred spatial distributions of surface crustal deformations are different for the seven slow-slip events. This phenomenon was also found by Kano et al. (2018) for five slow-slip events in the source area beneath Yaeyama Islands along the southwestern Ryuku subduction zone of Japan. This phenomenon might be due to different nucleation styles at different localities in the same source area. Nevertheless, more studies are necessary to understand the mechanism that generate slow-slip events, because they are still not known (e.g., Cook 2019). Xia et al. (2011) found infrasound signals recorded at Beijing about 3 days before the mainshock. This is the unique report of infrasound before the event. In fact, there are numerous infrasound recording stations around the world, especially at the volcanic areas in Japan which is somewhat close to Taiwan. Why were not the infrasound signals recorded at any station in Japan?

Infrasound
The effective infrasound speed can be written as (cf. Negraru et al. 2010): v i = v s + η·υ w , = v s + |υ w |cos(θ) where v s is the sound speed, υ w is the wind velocity with a magnitude of |υ w |, η is a unit vector in the direction of infrasound propagation, and θ is the intersection angle between η and υ w . Hence, v i is between v s -|υ w | and v s + |υ w | depending on θ. The sound speed, v s , in an ideal gas under adiabatic conditions, is given by: v s = ζRT K where ζ is the adiabatic index (= C p /C v where C p and C v are, respectively, the specific heat at constant pressure and that at constant volume), R is the gas constant, and T K is the absolute temperature. Clearly, the infrasound propagation is controlled by the three-dimensional spatial distributions of temperature and wind velocities in the atmosphere. Numerous studies (e.g., Drob et al. 2003;Matoza et al. 2017) show that the infrasound signals generated from a source are not spatially uniform and timevarying. Hence, it is reasonable that the infrasound caused by an earthquake may be detected at some places and not at others. Other factors, including natural events, animal communication, human-created noise, etc., can also yield the infrasound. A lack of both infrasound records at other stations and complete data of atmospheric physical conditions in eastern Asia before the earthquake makes us unable to confirm if the infrasound signals reported by Xia et al. (2011) is a reliable precursor or not.
In principal, pre-seismic slip could be a main factor in controlling the generation of infrasound signals. From the data shown by Xia et al. (2011) for global events, we cannot see a positive correlation between the precursor time and the earthquake magnitude. And the precursor times of some events are longer than 7 days. The two kinds of information might indicate that time-increasing preseismic slip plays a minor role on the generation of infrasound signals. Liu et al. (2001Liu et al. ( , 2004a suggested that TEC significantly decreased around the source area in the afternoons of the 4 days, 3 days, and 1 day before the mainshock (see Fig. 6a). Chuo et al. (2002) showed that the perturbation appeared 3 -4 days before the mainshock. Their observations also exhibit significant decreases in TEC and f o F 2 on several days in much earlier time before the mainshock. Heki and his co-authors (Heki 2011;Cahyadi and Heki 2013;Heki and Enomoto 2015) found that TEC anomalies started 40 to 70 minutes before several large earthquakes with M w ≥ 8.2 and reached nearly ten percent of the background TEC. They also stressed that the TEC amplitude depends on earthquake magnitude. Through rigorous examinations, He and Heki (2017) analyzed vertical total electron contents (VTEC) recorded near the epicenters before and after 32 earthquakes with M w 7.0 -8.0. They found that eight earthquakes with M w 7.3 -7.8 showed preseismic anomalies starting 10 -20 minutes before the events. They concluded that preseismic TEC anomalies occur only for few M w ≤ 8 earthquakes and their precursor times are shorter than those for M w > 8 events. The general pattern of temporal variation of TEC changes from their observations is shown in Fig. 6b. Since the coseismic anomalies are higher than the preseismic ones for most earthquakes, the peak value may appear at t > t r as displayed in Fig. 6b. Unlike Fig. 6a, the time scale of Fig. 6b is 'minute' rather than 'day.' In this figure the precursor time is T = t r -t o where t r and t o are, respectively, the occurrence time of an earthquake and the starting time of TEC anomalies. A comparison between Figs. 6b and 3 suggest that TEC anomalies appear only in the latest time span of step 3c. Figure 6a exhibits that the TEC anomalies immediately before and after the mainshock are very small. This is quite different from Fig. 6b. Wang (2021) proposed a model to link the electric field generated on the ground surface to the slip on a fault. When the slip of microfracturing on a fault becomes large enough to generate strong ground electric charges/currents, the EM precursors including the TEC anomalies may be produced. From Fig. 3, large enough slip can exist only in step 3c in a short time before a large earthquake. From core samples drilled from boreholes cutting through several faults on which large earthquake happened, Wang (2021) found that the thicknesses of fault core and gouges and the electric resistivity of gouges are important factors in controlling the generation of electric field on the ground surface. Among six cases, only the fault for the M w 9.0 Tohoku-Oki M w 9.0 earthquake of 11 March 2011 had the conditions to generate the ground electric field because its thicknesses of fault core and gouges are thick enough (about 10 m and 4.86 m for the former and the latter, respectively) and its resistivity is low enough (< 2.5 Ω-m). For the Chelungpu fault, the thickness of fault core, thickness of gouges, ad electric resistivity are 1.2 m, 1.12 m, and 6.9 Ω-m, respectively. Clearly, the conditions of the Chelungpu fault are not good enough to generate a strong ground electric field. In other word, the TEC anomalies either could not be generated from preseismic slip or were very weak before the 1999 Chi-Chi earthquake. Moreover, the TEC anomalies shown in Heki and his co-authors are always higher than the reference level; while those given by Liu's group are always lower than the reference one. It is difficult to explain this difference based on the positive correlation between ground electric field and slip obtained by Wang (2021). Note that the Chi-Chi earthquake was not included in the list of events having TEC anomalies through rigorous examination made by He and Heki (2017). Hence, it is necessary to deeply explore the difference in the pre-seismic TEC anomalies between Liu's group and Heki's. Tsai et al. (2006) found that anomalous lightning occurred 4 days before the mainshock. From Liu et al. (2001), the NmF 2 (= f o F 2 /80.3) that is the electron density at the F peak and TEC were anomalously smaller than the lower bound on 17 and 18 September 1999, which were 4 and 3 days before the mainshock, respectively. Hence, Tsai et al. (2006) assumed that the coincidental occurrence in time of enhanced cloud-to-ground lightening and reduced ionospheric TEC appears to suggest a possible coupling mechanism between these two phenomena. But, this assumption is questionable. On 16 and 19 September 1999, both lightning and TEC were strong and approached the upper bound. On 20 September 1999 which was just one day before the mainshock, there was no lightning but TEC appeared even though it was weak and closed to the lower bound. There is inconsistence between the observations in two time intervals. This implicates that the assumption given by Tsai et al. (2006) is questionable. Ohta et al. (2001) observed the changes in ELF emissions at Japan about 4 -5 hours before the mainshock. This anomaly observed in Japan appeared very later than those observed by Akinaga et al. (2001) done in Taiwan. Hence, it is necessary to explore the reason why there is such a big time difference for a similar EM phenomenon. Hayakawa (2013) stressed that EM emissions may abruptly increase a few hours before the earthquake. The observation in Japan by Ohta et al. (2001) could be the 'abrupt increase' in EM emissions. Hence, their observation could be reasonable. Nevertheless, a question arises: Why did not Akinaga et al. (2001) report such an 'abrupt increase' from the records in Taiwan? A possible reason is that the frequency bands are different between the two studies: 3 -30 Hz for Ohta et al. (2001) and 300 -3000 Hz for Akinaga et al. (2001). The former is much lower than the latter. A problem arises: Is it easier to generate lower-frequency EM emissions than to generate higher-frequency ones from the fault zone in a few hours before an earthquake? EM emissions are generated by ground electric fields caused by slip of microfracturing in a fault zone. From Fig. 3, the slip u due to microfracturing in step 3c for producing ELF signatures observed by Ohta et al. (2001) is longer than that in the earlier time of step 3b for yielding ULF signatures done by Akinaga et al. (2001). This might suggest that longer u leads to longerwavelength or lower-frequency EM emissions; while shorter u results in shorter-wavelength or higher-frequency EM emissions. In other words, higher-frequency EM emissions transferred to lower-frequency EM emissions in a few hours before the occurrence of an earthquake. This suggests that step 3c started a few hours before an earthquake. This may also answer the question why the UFL emissions reported by Akinaga et al. (2001) disappeared a few hours before the Chi-Chi earthquake.

ELF/ULF Emissions
Nevertheless, two problems should be resolved in the future. The first problem is: Why did the anomalous ELF signatures reported by Ohta et al. (2001) appear only in a short time span about 4 hours and then disappeared before the earthquake? The anomalous ELF signatures should continue until the earthquake occurrence because the preseismic slip due to microfracturing in the Chelungpu fault increased with time until the failure. The EFL emission should be recorded during the earthquake rupture because co-seismic slip is much longer than pre-seismic slip. The second problem is: Why were the anomalous ELF signatures not detected at other stations? Based on the model proposed by Wang (2021), a complete data set including different EM signals and TEC anomalies will help earthquake scientists to infer the accelerating process of slip in few days or a few hours before an impending earthquake.

Earthquake Lights
Earthquake light may occur before and during an earthquake (Derr 1973). Several mechanisms, including frictional heating to incandescence, conductivity of rocks, sparks caused by the piezoelectric effect, streaming potential, charge separation, bombardment between plasmas and charged particle, etc., were proposed to interpret the generation of earthquake light (see Derr 1973Derr , 1986Lockner et al. 1983;Lockner and Byerlee 1985). Lockner et al. (1983) and Lockner and Byerlee (1985) preferred to the piezoelectric effect. They showed that rock fracturing under water produces spectra of both atomic and molecular hydrogen, and from this and other reasons they concluded that earthquake lights are caused by week electron excitation of the ambient medium. Based on the model proposed by Wang (2021), earthquake lights must be stronger during the mainshock than before it, because co-seismic slip was longer than pre-seismic one. The frequency of observed earthquake lights is usually higher during and after an earthquake than before it (Derr 1973). However, Chen et al. (2000) reported that the co-seismic earthquake light was seen only once by local people at a site in the northern segment of the fault. This might be due to a possibility that during the earthquake most of the local people living near the fault only took care of themselves and thus ignored the nearby natural phenomenon. Song et al. (2006) observed that the concentrations of chloride abruptly increased on 19 September about two days prior to the mainshock and lasted a few days afterward at both hot and artesian springs at Kuantzeling, Chiayi (see Fig. 1) in west-central Taiwan. The temporal variation in spatial distributions of daily strain rates inferred by Chen et al. (2015) shows that the strain rate in the area, including Kuantzeling, to the southeast of the Chelungpu fault had become negative since 50 days before the mainshock. This means that the loading stresses on the underground rocks underneath Kuantzeling had been compressive since 50 days before the mainshock and lasted until the earthquake happened. The compressive stresses might decrease the porosity of affected strata, thus pushing the fluids in cracks flow into the aquifer. Upwelling of groundwater from the deeper rocks might carry the chemical components originally existed at depths up to the shallow part, thus leading to the mixing of deeper aquifer and causing groundwater chemistry anomalies. This leads to the anomalous increases in chloride concentrations on 19 September at Kuantzeling.

Groundwater Chemical Anomalies
A question arises: Is this anomaly a chemical precursor of the 1999 Chi-Chi earthquake? Unlike the Puli area, Kuantzeling is > 60 km to the south of the mainshock epicenter and outside the deformed zone of the fault. The preseismic effect should be much lower at Kuantzeling than at Puli. The types of anions with anomalous concentrations were different in the two areas: SO 4 2and NO 3at Puli and Clat Kuantzeling. In addition, the anomalies appeared about 200 days before the mainshock in the Puli area, yet only 2 days at Kuantzeling. This time difference is very big. I assume that the anomalous concentration of chloride appeared on 19 September might be associated with other events occurring in the Chiayi and Tainan areas rather than with the Chi-Chi earthquake. Chen et al. (2000) found abnormal behavior for 12 kinds of animals about few hours to 10 weeks before the mainshock (see Table 2). This means that such anomalous animal activities started from the later time of step 3a of Stage 3 (see Fig. 3). It is questionable that ants built new nests on trees about 8 -10 weeks before the mainshock. The anomalous activities of ants occurred were 2 -3 months before the mainshock, in other word, in the summer time from May to July. In summer, typhoons and storms usually lead to heavy rain in Taiwan. Meanwhile, the locality of this report is to the north of the Chelungpu fault and near the Dajia River. So, the weather condition may play a more important factor on such anomalous activities of ants than pre-seismic slip. Hence, I wonder if this abnormal phenomenon is an earthquake precursor or not. Of course, Table 2 exhibits that the anomalous activities of other animals occurred in an acceptable time period. Hence, it is significant to explore the possible reasons to cause those activities.

Biological Precursors: Animal Anomalies
Pre-seismic animal anomalies have been reported for a long time in both historical documents and scientific reports. Nevertheless, it is necessary to examine the reliability of observations and to explore the possible mechanisms. Several authors (e.g., Buskirk et al. 1981;Bhargava et al. 2009;Fan 2018;Woith et al. 2018) reviewed anomalous animal activities before large earthquakes. Buskirk et al. (1981) reported in details the anomalous animals in different epicentral distances and time intervals. Although there are numerous possible reasons to make animals living in the source area became abnormal before a large earthquake, the real reasons are not yet completely known. Hayakawa (2013) summarized the possible factors in causing abnormal animal activities from numerous studies and reports. They are: (1) a change in atmospheric pressure, (2) a change in gravity, (3) ground deformations (including uplift and tilt change), (4) acoustic signals and vibrations due to microcracks, (5) electromagnetic effects, (6) a groundwater level change, and (7) emanation of gases and chemical substances. Rikitake (1998) preferred the electromagnetic effects, though he also considered the importance of fourth and seventh factors. Grant et al. (2011) took the seventh factor to be the major one. Fidani et al. (2014) and Grant et al. (2015) considered local air ionization (belonging to the fifth factor) caused by stress-activated positive holes is one of the main reasons, especially for cows. Rikitake (1998) first related the epicentral distance, D, of an anomalous animal activity before an earthquake to its magnitude, M in the following relationship: .
. ( ) log M D 1 86 2 6 = + From the plot in Rikitake (1998), D is almost the upper bound epicentral distance. For the 1999 M 7.6 Chi-Chi earthquake, Eq. (5) leads to D = 161 km. This means that the animals living in an epicentral distance ≤ 161 km could become anomalous before the event. He also plotted the data points of T versus M. However, the value of T is distributed over a wide range from a few minutes to hundreds of days for each M. Hence, he did not inferred the relationship between T and M. From the occurrence histogram of log(T) (in units of days), he suggested that T is mainly in a range 1 -10 days. For the Chi-Chi earthquake, except for ants the values of T for other animals (see Table 2) are within this range and thus acceptable. Hayakawa (2013) utilized three plots to characterize the unusual animal behavior: (1) earthquake magnitude (M) versus D; (2) log(T) versus M; and (3) occurrence histogram of log(T). He assumed that these plots are comparable with the corresponding plots for different seismo-electromagnetic effects (radio emissions in different frequency ranges, seismo-atmospheric and -ionospheric perturbations). He also claimed that ULF and ELF electromagnetic emissions exhibit a very similar temporal evolution with that of anomalous animal behavior. This means that ULF and ELF electromagnetic emissions could play a primary role on inducing anomalous animal behavior. It is also suggested that a quantity of field intensity multiplied by the persistent time (or duration) of noise would significantly produce anomalous animal activities before an earthquake. Table 1 shows that both geomagnetic anomalies and anomalous ULF signatures appeared about 2 months before the mainshock and lasted until its occurrence. This time interval is comparable with those of anomalous animal activities as listed in Table 2. Hence, either geomagnetic anomalies or ULF signatures might be a significant factor in producing anomalous animal activities.

Some Pre-Seismic Natural Phenomena
Among non-animal precursors, sky luminescence was observed before several earthquakes, such as the 1965 Matsushiro swarm in Japan, the 1973 Veracruz event in Mexico, and the 1977 destructive Vrancea event in Romania (Derr 1973;Lomnitz 1994). Thunder-like sounds, changes in the level of ground water, and an increase in radon were reported as either coseismic phenomena or precursors of the 1975 Haicheng earthquake and the 1976 Sonpan earthquake (Deng et al. 1981;Teng and Henyey 1981).  found the occurrences of strong winds, skylight, thunder-like sounds and unusual gaseous odors before the mainshock. They assumed that these signals might equally serve as indications of an imminent earthquake as those based on aberrant animal behavior. But, some of them are not so reliable and more studies are necessary to examine the observations.

FUTURE CHALLENGES
This study reveals that most of observed precursors are reliable even though few of them must be re-examined. This suggests the possibility of earthquake prediction. Although reliable precursors may give us a clue to judge if an earthquake will happen in an area or not, observed data themselves cannot directly predict anything. We may learn something from a famous historical example. From a large number of high-quality astronomical data observed almost in sixty years by an excellent Danish group led by Tycho Brahe (1546 -1601), Germany astrologist Johannes Kepler (German, 1571 -1630) reduced three empirical laws of planetary motions. This is a great work providing inspirational experience. However, astrologists cannot predict the location of a planet at a certain time just based on the three empirical laws. Later, from Kepler's laws British physicist Issac Newton (1642 -1727) proposed the gravitational law that is a physical model based on rigorous mathematics. Astrologists can predict the location of a planet at a certain time from Newton's law. This example reveals that scientific revolution is like the motion of a science bicycle having two strong wheels: one for reduction (leading to empirical laws from observations) and the other for deduction (resulting in mathematical models). New observations may examine the existed models; while new models may test observations. When one of the two wheels is much larger than the other, the bicycle cannot move well. (Of course, this might be fine for a very skillful bicycle-driver.) We may learn something significant from Kuhn (1962) who wrote an excellent book to explain the structure of scientific revolutions.
For the earthquake prediction research, the wheel for observations (reduction) is much larger than that for theories (deduction). This cannot make earthquake prediction be successful. Hence, earthquake scientists should adjust the two wheels. One of the most important things is that earthquake scientists must try to construct the physical plus chemical models for respective precursors or even a unified model for all precursors. Based on the models or the unified model, earthquake scientists may predict an earthquake. In other word, earthquake scientists still have to face numerous strict challenges in the front of them.
In addition, numerous observations (e.g., Wang 2018a), laboratory experiments (e.g., Dieterich 1979;Ruina 1983;Marone 1998), and theoretical studies (e.g., Keilis-Borok 2002;Rice 2006;Wang 2009Wang , 2017 all suggest that faulting behaves like a non-linear process. This means that the earthquake rupture processes could be highly sensitive to the initial physical and chemical conditions of the crust and upper mantle in a large range around the hypocenter before an earthquake. This decreases the possibility of successfully predicting an earthquake and increases the difficulty of constructing an acceptable prediction model. One of the significant ways is to re-do the simulations when new observations are made.
Meteorologists have developed acceptable models for local and global weather forecasting, even though each major country has her own models. Although up to date weather forecasting is not 100% successful, it is still very useful for human activities and able to reduce hazards. Earthquake scientists must learn something from meteorologists and construct acceptable prediction models in the near future.

CONCLUSIONS
Earthquake prediction has been a long-term debatable problem in earthquake science. In order to resolve the problem, one of the ways is to study the possible precursors of a single earthquake. Numerous precursors were studied after the M s 7.6 Chi-Chi, Taiwan, earthquake of 20 September 1999. This makes us a chance to explore the debatable problem. Based on the time window (or the precursor time, T), earthquake prediction is classified as: very-long-term prediction (T > several ten or hundred years); long-term prediction (ten years > T > three years); intermediate-term prediction (T = six months to three years); short-term prediction (T = eight days to six months); and imminent prediction (T ≤ seven days). Meanwhile, all given precursors are classified into four categories: mechanical, electromagnetic, geochemical, and biological precursors. Each category may include several items. All given precursors with their respective precursor times are compiled. These precursors have been examined based on known physical and chemical theories. Moreover, the possible correlations between two precursors or among several precursors have also been discussed. Results reveal that most of observed precursors are reliable even though few of them must be re-examined. This suggests the possibility of predicting an earthquake. Nevertheless, up to date earthquake scientists cannot predict an earthquake just based on observations. Like weather forecasting, it is very necessary to construct the physical and chemical models of respective precursors or to develop a unified model for all precursors. Based on the models or a unified model, earthquake scientists may predict an earthquake. This means that earthquake scientists still have to face numerous challenges in the future.