Correlation Between Time Change in Modulus of Short-Period Geomagnetic Variation and Seismicity in Taiwan

1 Department of Earth Sciences, National Taiwan Normal University, Taipei, Taiwan, ROC 2 Department of Earth Sciences, National Chengkung University, Tainan, Taiwan, ROC 3 Institute of Earth Sciences, Academia Sinica, Taipei, Taiwan, ROC * Corresponding author address: Prof. Kuang-Jung Chen, Department of Earth Sciences, National Taiwan Normal University, Taipei, Taiwan, ROC; E-mail: kjchen@ntnu.edu.tw doi: 10.3319/TAO.2007.18.3.577(T) In this study, geomagnetic data of the Lunping observatory from 1993 to 2000 are utilized for computing the amplitude variation of short-period geomagnetic total intensity data, using the complex demodulation method (CD method). In order to compare these time changes with seismicity, earthquakes that occurred within 150 km of Lunping, with magnitude ML greater than 3.0, are located. The total sum of those earthquakes, summed month by month, is correlated with the modulus. After removing seasonal effect, our results show that the modulus of periods 24, 12, and 8 hr reveals a notable change that seems to be related to the total sum of events within the whole study period. One possible precursor is found 6 months prior to the 1999 high seismicity. The modulus for the periods 24, 12, and 8 hr increased gradually from the beginning of 1999 to August 1999. After earthquake occurrence the modulus decreased again to a normal level. We propose that this notable increase might be related to a preparation process for this strong earthquake. (


INTRODUCTION
first suggested the detection of stress-induced changes in local geomag-netic fields as a means of monitoring tectonic activity.Laboratory work accomplished by Kalashnikov (1954), Kapitsa (1955), and Ohnaka and Kinoshita (1968) has shown that a change in stress can produce changes in magnetic field.The theory of stress-induced changes in rock magnetization is generally in agreement with laboratory studies (Stacey 1964;Nagata 1970;Stacey and Johnston 1972).Anomalous electromagnetic precursors to earthquake have been reported in many papers (e.g., Beamish 1982;Rikitake 1987a, b;Fraser-Smith et al. 1990;Yoshino et al. 1993;Oike and Murakami 1993;Fujinawa and Takahashi 1994;Pulinets et al. 1998;Rikitake 1997;Liu et al. 2000;Zeng et al. 2001;Nagao et al. 2002;Zeng et al. 2002;Yen et al. 2004).Research on the relation between time change of conductivity and seismic activity is an important subject in earthquake prediction.One well-known model of an earthquake mechanism is called the 'dilatancy' model, which proposes that volume increase in stressed rock supposes a large increase in porosity and mechanical permeability (Scholtz and Kranz 1974).When electrolytic solution permeates the pores and cracks, the conductivity of rock can increase significantly.Apart from the 'dilatancy' model, there are other reasons, such as a rise in temperature, which may also lead to noticeable increase in conductivity within the rock in an earthquake preparation zone (Olhoeft 1981).Jones and Price (1970), and Bailey and Edwards (1976) developed a model to simulate geomagnetic variations on the Earth's surface caused by underground anomalous conductivity.Hautot (2005) developed a model to simulate geomagnetic variations on the Earth's surface caused by underground anomalous conductivity.Electromagnetic phenomena are now considered a promising candidate for short-term prediction of high seismicity (e.g., Hayakawa and Fujinawa 1994;Hayakawa 1999;Hayakawa and Molchanov 2002).
The CD method is one of the Fourier-based time series analyses.Some tedious computation in time domain can be simplified by the CD method, which performs the operations in the frequency domain.Banks (1975) first introduced the application of the CD method into the field of geophysics.Chen et al. (1985) applied a similar CD method to geomagnetic data in Taiwan.The present paper, which describes the CD technique and its computational implementation, relies heavily on that work.In this study, we have designed a computer program for the CD method based on formulations by Banks (1975).Because the Taiwan region is an active seismic zone, the geomagnetic data from 1993 -2000 of Lunping Observatory (121°10'E, 25°00'N, in Taiwan, see Fig. 1) were used in order to get a possible relation between the time change in modulus of short-period magnetic variations and seismic energy released to search for a possible precursor to earthquake occurrences.The instrument installed at this station is a continuous-recording proton precession magnetometer having an absolute sensitivity of 1.0 gamma (Geometrics Model G856).The proton precession magnetometer differs from classical magnetometers in that it relies on an atomic constant, i.e., the magnetic moment of a proton or hydrogen nucleus.Accordingly, a proton precession magnetometer is free from instrument drift, and is not affected by environmental factors such as temperature and humidity.The precision of the observation of the total intensity of magnetic field is within ±1.0 gamma ( γ ).The sensor of the magnetometer is supported by a wood-plugged PVC tube and protected from rainfall and dust by a larger PVC tube (Huang 1990).Sampling is made at an interval of 1-minute.

The Definition of Complex Demodulation
The process of complex demodulation, as applied to a time series X(t) (assumed to be sampled at a series of points spaced at equal intervals of time ( ∆t )), is defined in terms of two very simple mathematical operations.First, each frequency band of interest in the spectrum is shifted to zero frequency by multiplying each term of the time series by a complex exponential function exp(-'t) iω , where ω' is the central frequency of the shifted band.A new series X S (ω',t) is produced for each frequency band, i.e.,: Second, the frequency-shifted series is then low-pass filtered using a set of weights a k (k = -m to +m, m = positive integer), and the result is the (complex) demodulated time series X ( ',t) d ω : X ( ',t) = a X ( ',t + k t) (2) The demodulates can be expressed most conveniently in the form: in terms of the modulus and phase of X ( ',t) d ω .

Demodulation of Data Containing Periodic Components
Let us suppose that the data contain a periodic component that would produce a peak in its Fourier spectrum at frequency ω o , i.e.,: (4) The simple complex demodulate centered on frequency ω ω δω ' = + will be: which contains frequencies −δω and − + ( ) 2ω δω o The frequency-shifted series is then low-pass filtered, so hopefully the component at − + ( ) 2ω δω is removed completely, leaving: Here we have assumed that the filter has unity response and introduces no phase shift at frequency −δω .
If we know that the data contain a periodic variation with a frequency of ω o , it is natural to choose ω ω ' = o as the central frequency of the demodulation, and to shift down to zero the band of frequencies immediately around ω o .In that case, δω = 0, and If A and γ change with time, there will be a corresponding change in the modulus and phase of the demodulation.

RESULTS
In order to identify the predominant frequency of short-period variation of geomagnetic total intensity, spectral analyses have previously been applied to the continuous records.Almost all the spectra in Lunping reveal that there are three relatively significant periods at 24, 12, and 8 hr.Two examples of the amplitude spectrum in Lunping station are shown in Fig. 2. Therefore, we choose the frequencies ω' = 1, 2, and 3 cycles/day (period = 24, 12, and 8 hr, respectively) as the central frequency for demodulation in this study.The CD method requires that the time series be continuous over a chosen time window.The results achieved by putting the whole time series into the CD method program are approximately similar to those obtained by dividing the time series into a number of segments, analyzing each segment and combining them together.Accordingly, yearly data sets from 1 January 1993 to 31 December 2000 are chosen to be processed by the CD method program.Figure 3 illustrates the time changes of modulus from 1993 to 2000 in periods 24, 12, and 8 hr, respectively.Gray dotted-lines represent the mean of each curve and its standard deviation.The bold dotted-line in each plot represents the linear regression.The linear regression uses the least squares method to construct a fit to a set of data points (x i , y i ) i = 1,...., n by a polynomial of degree p where: For comparison, we introduce the time changes in modulus observed at ASP observatory (Alice Springs, 133°53'E, 23°46'S), in a "seismic quiet" area (earthquakes occurred with magnitude greater than M L 3.0, within 200 km from these stations, are listed in Table 1), analyzed by the same procedure as mentioned above.The result is shown in Fig. 4. In general, the time changes of modulus contain a significant period of one-year, which can also be found in LNP.It seems that no abrupt increase is shown before the 1999 high seismicity, as the Fig. 3 shows.This indicates that the abrupt increase in LNP in unique.
In order to search for the correlation between the modulus of short-period magnetic variation and seismicity, earthquakes of magnitudes greater than 3.0 within a radius of 250 km of the Lunping Magnetic Observatory that occurred during the study period were collected from the Earthquake Catalogue of Central Weather Bureau.Total sums of earthquakes with magnitude (a) M L 3.0 and (b) M L 4.0, within a radius of 50, 100, 150, 200, and 250 km, respectively (top to bottom), from the Lunping Observatory were performed (Fig. 5).The hypocenter distance of the 1999 Chi-Chi earthquake for Lunping is about 160 km.This implies that the seismicity affected by the 1999 Chi-Chi earthquake can be neglected if we set 150 km as the criterion for choosing the earthquake.
The time changes of modulus shown in Fig. 3 seem to be ascribed to seasonal change.The seasonal effect can also be found in many areas (Shiraki 1980;Chen 1981;Gong 1985;Chen and Fung 1993).In order to remove the seasonal influence, the values of the same day were averaged over 8 years from 1993 to 2000 to get the annually variation.Figure 6 shows the result after the removal of seasonal effect.Positive anomalous variation (anomaly 1) can still be found appearing before the 1999 high seismicity.The amount of deviation in period 12 hr is about 2 and 1.5 nT in period 8 hr.Since the standard deviations of these two periods are 0.7746 and 0.5049 nT, respectively, the variation amounts above are larger than three standard deviations.However, after the earthquake, the variations are reduced to normal levels.Importantly, prior to the 1999 high seismicity, the regressions of modulus (solid dots in each curve) in period 12 and 8 hr show a significant increase (as arrows shown in Fig. 6), with about 3 nT, and it is reduced to its normal level later.We consider these variations as occurring locally and may be related to the 1999 high seismicity.

DISCUSSIONS AND CONCLUSIONS
Chen et al. ( 2006) found that time variations in transfer functions are generally correlated to seismic energy locally released, hence they concluded that transfer functions are adequate for short-term earthquake prediction.They also pointed out that the transfer functions of shortperiod might be meaningful for monitoring earthquake occurrence.In this study, analyzing the short-period of geomagnetic data, it was discovered that three significant frequencies (1, 2, and 3 cycles/day) in Lunping observatory are revealed.Although the detailed mechanisms are not fully understood, the authors speculated that they must be due to some source fields under the ground, somewhere around Lunping.However, more observations should be accumulated in order to confirm this phenomenon.
The time changes of the modulus of these three characteristic periods, observed at the Lunping observatory, show notable anomalous variation before the 1999 high seismicity and are reduced to normal levels after then.After removing seasonal effect and comparing results to those of a remote station (ASP)(a seismically quiet area), we can definitely conclude that a notable anomalous variation (Anomaly 1 in Fig. 6) appeared about 6 months before the 1999 high seismicity (S1 in Fig. 6).A recovery trend occurs after then.The global geomagnetic activity index (Kp), (acquired from National Geophysical Data Center (NGDC) for the study period), is plotted in the bottom of Fig. 6.It can be easily found that there is no notable correlation between Kp and the seismicity in Lunping area.It can be concluded that the variation in Lunping is a local effect rather than a global effect.
The selection of nearby earthquakes is an important matter, because so many earthquakes occurred in the Taiwan region during the study period (as Fig. 1 shown).Chen and Fung (1985) concluded that imaginary Parkinson arrows can be affected by geomagnetic induction within 150 km.The selection principle used here was to choose those earthquakes as near as possible (not farther than 150 km) and as large as possible, but this is certainly not objective enough.The difference in seismicity between that of magnitude larger than 3.0 and larger than 4.0 of the earthquake catalogue (released by CWB) within the study period is quite small (as Fig. 5 shown).If the distance is extended to 250 km from Lunping, the conclusion that a correlation exists between modulus and seismicity is still valid.

Fig. 1 .
Fig. 1.The location of the Lunping Observatory and the earthquakes occurring within the interval 1992 -2001, with a magnitude M L greater than 3.0 within 150 km distance from Lunping.Light gray open circles denote the events with magnitude M L greater than 3.5 occurring in the Taiwan area in the same period.The epicenter of the 1999 Chi-Chi earthquake and the high seismicity are also shown.

Fig. 3 .
Fig. 3. Time changes in modulus of (a) period = 24 hr (b) period = 12 hr, and (c) period = 8 hr for the Lunping observatory in the interval 1993 -2000.The dotted lines represent mean, and long-dash lines denote its linear regressions.The arrow represents the 1999 high seismicity.
Table1.The earthquakes occurring within 200 km from stations ASP, during 1992 to 2001 with magnitude M L 3.0 (located by the Incorporated Research Institutions for Seismology (IRIS)).

Fig. 4 .
Fig. 4. Time changes in modulus of (a) period = 24 hr, (b) period = 12 hr, and (c) period = 8 hr for the Alice Springs' observatory in the interval 1993 -2000.The dotted lines represent mean, and long-dash lines denote its linear regressions.The arrow represents the 1999 high seismicity.

Fig. 5 .
Fig. 5. Time changes of monthly summation of earthquakes of a magnitude (M L ) (a) greater than 3.0 (b) greater than 4.0, within a radius of 50, 100, 150, 200, and 250 km, respectively (from top to bottom), for the Lunping Observatory.