The Characteristics of Magnetic Signals Observed at Lulin ELF Station and Their Association with Earthquakes in Taiwan

Electromagnetic wave data from years 2009 2010 which were recorded by the extremely-low-frequency (ELF) station at Lulin Observatory in Taiwan (120.87°E, 23.46°N; 2862 m) are analyzed to explore the wave properties and their association with seismic and geomagnetic activities. Specifically, the power spectra, the ellipticity and the orientation of the below-10-Hz wave signals are studied. The correlation coefficients of these ELF signals and geomagnetic activity indices Dst and Kp are within -0.2 to 0.2, indicating geomagnetic activities are not likely sources for these emissions. In the frequency versus time (f-t) spectrograms, a few events exhibit a frequency shift between 4 and 6 Hz and the possible association with seismic activities of magnitudes (M) larger than 4 in the vicinity of Taiwan are examined. The subset of earthquakes in Taiwan with M/r2 (r is the epicenter distance in kilometer from the Lulin station) of greater than 6 × 10-4, a total number of 50, show a higher association with ELF wave events showing a frequency down-shift within 7 days prior to the earthquakes and a frequency bounce-back within 7 days after the earthquakes. The epicenters of these earthquakes are all within a distance of 100 km from the Lulin station and with a depth of 50 km, suggesting that the detected ELF frequency-shift wave events maybe are pre-earthquake electromagnetic signals.


InTrOduCTIOn
Over the past decades, there have been many observation reports of electromagnetic disturbances in the frequency range of 10 -3 to 10 8 Hz before large earthquakes (e.g., Kalashnikov 1954;Johnston 1989;Parrot 1990;Hayakawa and Molchanov 2008).Among them, disturbance signals at the ELF range below 10 Hz were featured in several studies.Using ground instruments measuring electromagnetic signals between 0.1 and 10 Hz, Belov et al. (1974) found unusual electromagnetic disturbances at 10 seconds before occurring of a M = 6 earthquake in Kamchatka located within 70 km of the instruments.Gokhberg et al. (1982) and Warwick et al. (1982) investigated the association between seismic activities and the recorded magnetic disturbances; they proposed that the compressed rocks may emit electromagnetic waves hours before the occurrence of earthquakes; however the pre-quake electromagnetic disturbances may not be universal.Fraser-Smith et al. (1990) found electromagnetic disturbances 3 hours before the Loma Prieta earthquake (M = 7.1) with an instrument that measured electromagnetic signals from 0.01 to 10 Hz.Adams (1990) noted that the intensity of measured E&M (electromagnetic) signals below 1 kHz increased 10 days prior to the three earthquakes occurred in the same month in California, including the Loma Prieta earthquake.Dea et al. (1991Dea et al. ( , 1993) ) used the E&M data between 0.1 and 10 Hz as well as 10 and 40 Hz to discover that that the intensity of the disturbed signals between 0.1 and 5 Hz increased one day prior to an earthquake (M = 4.6) which occurred 200 km away in the southern California.Molchanov et al. (1992) discovered unusual electromagnetic signals at a few hours before the earthquakes of Loma Prieta and Spitak (M = 6.9).Kopytenko et al. (1993) analyzed earthquakes in Spitak and Japan (M = 6.4) using observational data from several stations and found this type of disturbances were not originated from magnetosphere.Schekotov et al. (2007) specifically analyzed 4 to 6 Hz radio frequency data to avoid interferences from global geomagnetic activities (frequencies < 2 Hz) and the background Schumann resonance signals (frequencies > 7 Hz); they found this band emissions in the Kamchatka region started 5 days before an earthquake and persisted for 5 days after the quake for the period from 24 February 2003 to 6 April 2003.In their work, they found, for the probable quake-related radio signals, the spectral power ratio of the radio signals recorded along the north-south and the eastwest directions exhibited enhancement, while the standard deviation of the ellipse orientation angle and the ellipticity decreased and thus the polarization of waves became more linearized.Recently, Freund et al. (2009) reported that a laboratory experiment involving rocks subjected to large stress discovered that field-ionization of air molecules and corona discharge at the rock-air interface can be analyzed indicating electromagnetic signals may be emitted by stressed-rocks before earthquakes.
In this paper, we analyze the 4 to 6 Hz magnetic field signals recorded by an ELF station in Taiwan between October 2009 and September 2010 to investigate possible associations with seismic activities (M ≥ 4) in the vicinity of Taiwan.

MAgnETIC ELF STATIOn, ELF SIgnALS And ThE dErIvEd WAvE PArAMETErS
With the goal of monitoring global lightning activity, a magnetic ELF station located at the Lulin Observatory-Yushan National Park (23.47°N, 120.87°E; 2862 m altitude; magnetic inclination: 16.76°N) has been installed and operated by the NCKU ISUAL team in Taiwan since August 2003 (Wang et al. 2005;Huang et al. 2011).The system comprises a pair of EMI-BF4 magnetic coils that are sensitive in the radio frequency band between 0.3 and 500 Hz.The orientations of the coils are either parallel (H; northsouth) or perpendicular to (D; east-west) the geomagnetic field.The system is also equipped with a GPS clock and all the data were sampled continuously at 5 kHz.In the past, in addition to performing routine recording of lightning sferics, intriguing ELF-whistlers in the frequency range between 60 and 100 Hz had also been observed by this system (Wang et al. 2005(Wang et al. , 2011)).
In this study the three primary wave parameters to be constructed from the magnetic ELF signals, as suggested by Schekotov et al. (2007), are the power spectra, the ellipticity (tanβ) and the angle (θ) between the major axis of the polarization ellipse and the H-axis.The power spectra of the waves include P hh (N-S component), P dd (E-W component) and P hd (cross spectrum of H and D components), the power ratio P hh /P dd , as well as the total power P hh + P dd .The helicity of the wave polarization traditionally is represented by the angle β as defined in Eq. ( 1), where β > 0 for righthand polarized waves and β < 0 for the left-hand polarized waves.
sin Im The ellipticity is defined by tanβ and it is the ratio between the minor axis and the major axis of the polarization ellipse.The angle θ between the major axis of the polarization ellipse and the H-axis can also be calculated.Angle θ satisfies the following equation: The frequency-time (f-t) spectrograms for the total power, the power ratio, tanβ, and θ can then be plotted.Figure 1 shows an example of the derived f-t spectrograms for the wave data recorded in the month of December 2009.A specific feature which can be immediately perceived from the total power spectrum (Fig. 1a) is that there is a strong emission band at 4 to 6 Hz.Also the peak intensity can shift between 4 and 6 Hz.The dates of the frequency-shift are marked in the Fig. 1c that contains the power ratio spectrum.The two dark-pink rectangles denote the dates when the wave frequency drops from 6 to 4 Hz which are identified to be December 8 and 29, respectively.The yellow rectangle denotes the date when the wave frequency bounces back from 4 to 6 Hz which is December 15. Figure 1b shows that the calculated spectral values of tanβ are between -1 and 1.The dark-pink rectangles mark the dates when the values approach 0, where the spectra are denoted by light green colors (see the color bar to the right of this panel) and the frequency is near 4 Hz.There are intermediate values between 4 to 6 Hz and for this moment we track 4 Hz in our analysis.Figure 1d shows the spectrum for θ, in which white rectangles denote those with θ > 0 (red color) among the θ < 0 emission band (blue color) while the dark-pink arrows indicate those with θ < 0 (blue color) in the θ > 0 emission band.The angle θ was found to concentrate more at a certain value before large quakes (Schekotov et al. 2007), so here we tag the sharp changes associated with focused values of θ in this rather continuous band.

COrrELATIOnS WITh gEOMAgnETIC AC-TIvITIES And SEISMIC ACTIvITIES
To assess the effects of geomagnetic activities on this observed 4 -6 Hz band, correlations between geomagnetic activities and these recorded signals at the frequency band under 10 Hz are further examined.The geomagnetic indices Dst and Kp are adopted since they are applicable for our low-latitudinal ELF station.Figure 2 shows the calculated correlation coefficients for the data in the month shown in Fig. 1.The values generally fall between -0.2 and 0.2 for the entire under-10-Hz frequency band, meaning little correlations between these indices and the detected waves below 10 Hz.Correlation coefficients calculated for the data collected at other months also exhibit low values.Lag correlations are further explored by calculating the cross-correlation coefficients specifically for the 4 and 6 Hz signal components and the results are shown in Fig. 3. Figures 3a to d are the correlation coefficients for the 4-Hz-wave component, and Figs.3e to h are those for the 6-Hz-wave component.The calculated cross-coefficients are small between the wave parameters and both Dst and Kp; they also are insensitive to length of the lag days (up to seven) as shown in the figures.
On the other hand, for the responses of the 4 -6 Hz waves to seismic activities, Schekotov et al. (2007) reported the observed waves showed enhancement in signal power, the waves became linearized (tanβ = 0), and the values of θ made sharp transitions between 4 and 6 Hz for their measured magnetic signals, from about 5 days before the earthquakes.Since our measured signals have rather special features between 4 and 6 Hz (Fig. 1), comparisons with seismic activities around Taiwan are then implemented.First, the earthquakes with magnitude M larger than 4 which occurred between October 2009 and September 2010 are collected from the archived database of the Central Weather Bureau (CWB) in Taiwan.The total number of M > 4 quakes is 169 during this period.Figures 4a, b, and c depict the magnitude M, the epicenter distance to the Lulin station r, and M/r 2 versus to time during the study period (October 2009 -September 2010), respectively.In the top three panels of Fig. 4, only 50 earthquakes, with M/r 2 larger than 6 × 10 -4 (km -2 ), are shown.The reason to plot the earthquake data above this value is outlined below.
For easy comparison, temporal distributions of the 4 -6 Hz wave events for the same period are shown in Figs.4d to e. Figure 4d denotes the time when the event frequency bounces from 4 to 6 Hz (the solid diamonds) and the time the wave frequency shifts from 6 down to 4 Hz (the hallow diamonds).Figure 4e presents the times when tanβ equals 0. Figure 4f shows the times where θ > 0 or θ < 0. The selection criteria are the same as those depicted in Figs.1b  and d.
To examine the association between two seemly randomly distributed lists of events is not trivial.It is especially true for seismic activities containing a two-dimensional distribution of space and time.Since we have focused our study on a specific region, we can treat the seismic events as a set of randomly-distributed events in time.In order to analyze the correlations between seismic activities and the observed 4 -6 Hz transition wave events, one can calculate the association number N (Hsu and McPherron 2002).
The association number is a generalized correlation analysis to examine two randomly distributed event lists.
In this analysis 2h is the temporal window size of the earthquake event, denoted as b j , u is the time shift (delay) between event b j and an ELF-wave event a i .If a i is within the window of b j after a time shift of u, then a i counts as an associated event.N will be the total number of the associated events.Different time delay u can give different N and there could be a value of time delay u that produces a maximum value of N, which suggests a possible lag relation between two events.The calculated association numbers N between earthquakes and the wave spectra and polarization parameters for the data shown in Fig. 4 are presented in Figs. 5 -7.The scan of the time delay embraces ±4 hours to ±7 days, with a step size of 4 hours.Figure 5 shows the association N between earthquakes, 4 to 6 Hz frequency up-shift (Fig. 5a) and 6 to 4 Hz frequency down-shift (Fig. 5b) of 4 -6 Hz events for the set of data shown in Fig. 4d.The black line denotes the association number for the time delay u, the light green line shows the smoothed values of the association number.The horizontal axis is the time in units of days, day 0 is the day when an earthquake occurs.Evidently, the number N peaks at around u = -2 days at the top panel and this indicates wave frequency down-shift mostly occurs at two days before the earthquakes.Yet, after the occurrence of earthquakes, the association number drops sharply, suggesting that the wave frequency down-shift can likely be found prior to the occurrence of earthquakes.In addition, the total association number within seven days prior to the quakes (the area under the blue curve) is 26.17, more than a factor of 2 higher than the value of 10.17 within seven days after the earthquakes.This also indicates statistically that the probability for frequency down-shift events to occur is higher within 7 days before an earthquake than after.Meanwhile, the bottom panel shows that wave events exhibiting frequency up-shift (bounce back) mostly occur ~4 to 5 days after the earthquake.The total association number shows an opposite trend compared to that in the top panel, the number is 10.58 within seven days before the quake, while the number increases to 27 within seven days after; also, the ratio is a factor of more than 2, indicating an increase in trend for frequency up-shift events to occur after an earthquake.Note that only earthquakes with M/r 2 larger than 6 × 10 -4 were selected to compute the association number, because the feature involving the frequency down-shift at two days before the quake and the frequency up-shift at 4 -5 days after the quake is consistently present for earthquakes with M/r 2 larger than 6 × 10 -4 .For earthquakes with M/r 2 smaller than this value, even though the total numbers of earthquakes increase substantially but the above mentioned feature no longer hold.Thus only earthquakes with M/r 2 larger than 6 × 10 -4 are plotted in Figs.4a to c. Hence M/r 2 is a simple representation of the intensity of seismic events around the Lulin ELF station.
Figure 6 shows the association number N for earthquakes with waves having tanβ = 0 in Fig. 4e.Two peaks occur before the occurrence of earthquakes appearing near u = -2 days and -5 days, but the difference in N between pre-earthquakes and post-earthquakes are not as distinct as those in Fig. 5. Figure 7 shows the association number N computed for wave events with θ < 0 (Fig. 7a) and θ > 0 (Fig. 7b) for the data shown in Fig. 4f.The association number peaks at longer time delay before the earthquakes than that after the quakes, but then again the feature is not as distinct as that shown in Fig. 5.

dISCuSSIOn
ELF magnetic wave data collected by the ELF station at Lulin Observatory between October 2009 and September 2010 were analyzed to investigate the features of waves with frequency below 10 Hz.Wave frequency shift between 4 and 6 Hz were identified and their correlations with geomagnetic activities were examined and found to be insignificant (Fig. 3).
Correlations of these wave events with earthquakes (M > 4) in the Taiwan region were studied for ±7 days within the occurrence of earthquakes.Frequency down-shift (6 to 4 Hz) wave events occurred often before these earthquakes while the frequency up-shift (4 to 6 Hz) wave events often occurred after earthquakes (Fig. 5) for the earthquake dataset with M/r 2 larger 6 × 10 -4 (km -2 ) (Fig. 4c).There are gaps in the earthquake database between days 47 to 71, 112 It should be noted that frequency shifts occurred less often (Fig. 4d) in these gaps in which earthquakes separated by more than 14 days.Also as Fig. 4 shows, a cluster of earthquakes can occur on the same day (e.g., day 154) and they can result in a high association number for that day.However, even we bin the earthquake data by counting earthquakes occurring on the same day as one earthquake and the peak values of N drops down, yet the total association number with wave frequency-shift events still shows similar trend before and after earthquakes (Fig. 8).Therefore, the above statistical result appears to be significant.If these wave events are indeed pre-earthquake electromagnetic signals (Freund et al. 2009), the electromagnetic signals emitted by the stressed rocks at 4 Hz may have been excited before the quakes.In fact, among the 50 earth-  (a) N with frequency shift-down events.This total association number is higher before earthquakes occur (Day = 0), with a value of 16.67, and it drops down to 10.17 after earthquakes occur.(b) N with frequency shift-up events.This number is higher after earthquakes occur (Day = 0), with a value of 16.17, while it is 9.0 before the occurrence of the earthquakes.quakes in the database, 37 are associated with thrust faults, accounting for 74% of the reported events.This type of fault is formed when crust is compressed.In addition, the epicenter distances of these earthquakes are all within 100 km from the Lulin station (Fig. 4b) and the depths are all within 50 km, therefore it is theoretically feasible that this enhanced wave emission from the nearby seismic activities can be more easily picked up by the ELF station, as illustrated by the epicenter distribution shown in Fig. 9.
Regarding the correlations of earthquakes with tanβ = 0, θ > 0, and θ < 0 events, Figs. 6 and 7 also show higher association numbers before earthquakes than those after.However, the difference is not that distinct as it has been mentioned above in the discussion related to Fig. 5. Therefore, a firm conclusion cannot be made here since it is possible that the select criteria to identify these events may alter the calculated association numbers.However, one conclusion still can be drew basing on results shown in Figs. 5  and 8, the opposite trend of frequency shifts at pre-and-post earthquakes are worth noting and further investigations using a longer term ELF data are clearly needed.

Fig. 1 .
Fig. 1.Frequency-time spectrograms derived from measured magnetic field signals on December 2009.The window size of FFT is 2 seconds.Spectra under 10 Hz are plotted.(a) Total power P hh + P dd .(b) tanβ.(c) Power ratio (P hh /P dd ) between N-S and E-W components.(d) θ.Schumann resonance around 8 Hz would be more distinguished with wider window sizes.Some emission lines crossing all frequencies are due to background noise enhancement.

Fig. 2 .
Fig. 2. Calculated correlation coefficients of signal parameters and geomagnetic index Dst and Kp under 10 Hz (horizontal axis) for December 2009.(a) Correlation coefficients of total power (red line) and power ratio (blue line) vs. Dst.(b) Correlation coefficients of total power (red line) and power ratio (blue line) vs. Kp.(c) Correlation coefficients of θ (red line) and tanβ (blue line) vs. Dst.(d) Correlation coefficients of θ (red line) and tanβ (blue line) vs. Kp.

Fig. 3 .Fig. 4 .
Fig. 3. Calculated cross correlation coefficients of signal parameters and geomagnetic index Dst and Kp under 10 Hz (horizontal axis) for December 2009.The lag days are from -7 to 7 days.(a) and (b): correlation coefficients for total power and power ratio for 4 Hz component vs. Dst and Kp, respectively.(c) and (d): coefficients for θ and tanβ for 4 Hz component vs. Dst and Kp, respectively.(e) to (h), same as (a) to (d) but for 6 Hz component.

Fig. 5 .
Fig. 5. Association number N for data of earthquakes and 4 -6 Hz events based on Fig. 4d.(a) N for frequency-down events.(b) N for frequency-up events.Seven days before the earthquake occurrence and seven days after earthquake occurrence are counted.

Fig. 9 .
Fig. 9. Geographical distributions of epicenters (green dots) of earthquakes with M/r 2 > 6 × 10 -4 (km -2 ) with respect to the Lulin station (blue triangle), the size of dots represents the magnitude of each earthquake as denoted at the top of this figure.