Azimuthal propagation of seismo-magnetic signals from large earthquakes in Taiwan

This study uses a network that is comprised of 10 total intensity magnetometers to detect azimuthal propagation of seismo-magnetic emission waves during 26 earthquakes that occurred between July 2007 and December 2008 in Taiwan. The propagation azimuth and phase velocity of the seismo-magnetic waves are calculated using frequency wavenumber analysis at the ultra low frequency of 0.05 Hz every 30 min. We superimpose the derived azimuths within a moving window of 30 days as the monitored distributions, and the entire dataset as the background distribution. We also find the propagation azimuths of the seismomagnetic anomalies of each earthquake by subtracting the background from the monitored distributions. The results show that frequency wavenumber analysis can be applied to evaluate azimuthal propagation of seismo-magnetic emission waves using a scalar of geomagnetic total intensity fields. The success detection rate of seismo-magnetic anomalies increases from 62% of the 26 earthquakes to 77% using the surface magnetic anomalous reference tip (SMART) to substitute the epicenters. Meanwhile, the odds proportions between the azimuths of the seismomagnetic emission waves towards and away from SMART reveal the associated anomalous propagation.


Introduction
Geomagnetic anomalies within a wide frequency range from direct current (DC) to very low frequency (VLF) that are associated with large earthquakes have been studied intensively [Hayakawa and Fujinawa 1994, Hayakawa 1999, Hayakawa and Molchanov 2002. Fraser-Smith et al. [1990] observed that amplitudes in the geomagnetic field in a frequency band between 0.05 Hz and 0.2 Hz suddenly increased a few hours before the Ms 7.1 Loma Prieta earthquake. Following their findings, geomagnetic anomalies in the ultra low frequency (ULF) range (#300 Hz) that are associated with large earthquakes have been examined [Bernardi et al. 1991, Molchanov et al. 1992, 1995, Kopytenko et al. 1993, Merzer and Klemperer 1997, Kawate et al. 1998, Gotoh et al. 2002, 2004a, b, Hattori 2004, Chen et al. 2011]. To locate a forthcoming earthquake, three-component geomagnetic data are analyzed using principal component analysis and signal value decomposition, to evaluate the major azimuth of the seismo-magnetic emission waves (SMEWs) [Gotoh et al. 2002, Hattori 2004, Telesca et al. 2004, Hattori et al. 2006]. Chen et al. [2009a] examined such magnetic anomalies statistically before 181 earthquakes (M w ≥4.3 or M L ≥5.0) and proposed the existence of seismo-magnetic fields. However, SMEW anomalies have not yet been tested in detail through any wave propagation analyses.
Seismometers with high sampling rates can precisely record changes in seismic waves. For teleseismic earthquakes, the delay time of seismic waves between every two stations can be accurately obtained by comparing the seismograms. To compute propagation directions and velocities of seismic waves, Capon et al. [1967] and Capon [1969a, b] developed frequency wavenumber (FK) analysis. This evaluates the directions and phase velocities from filtered seismic data, where the distance between every two stations is smaller than the wavelength of a single wave (i.e. a wave phase from 0 to 2r). Therefore, if SMEWs are indeed generated by an earthquake during a seismogenic process, it would be possible to estimate the direction of a forthcoming earthquake using the FK analysis.
In the present study, standard FK analysis is applied to analyze geomagnetic total intensity fields, which are scalars recorded by a network of 10 magnetometers (Table 1) in Taiwan, to determine the directions of SMEWs before large earthquakes (Figure 1). A statistical study of 26 earthquakes (M w >4.3; Table 2) is also carried out, which were retrieved from the broadband array in Taiwan , Liang et al. 2003, 2004 for the seismology for the period from July 2007 to December 2008.

Theory and the FK method
It has been known for many years that electromagnetic waves propagate at the speed of light of about 3 ×10 5 km/s in the air [Cheng 1989]. Since the phase speed is a product of frequency and wavelength, the wavelength of anomalous waves at 0.05 Hz will be about 6 ×10 6 km (as 3 ×10 5 / 0.05).
64 CHEN ET AL.  A phase change from 0 to 2r of such a signal wave will take approximately 20 s. Therefore, if the sampling interval of a magnetometer is shorter than 10 s, seismo-magnetic anomalies at about 0.05 Hz will be detectable. As spherical waves are transferred into plane waves due to the long propagation distance from teleseismic earthquakes, FK analysis can be applied to study seismic waves [Goldstein and Archuleta 1991a, b, Satoh et al. 2001, Schisselé et al. 2004]. The maximum power spectrum, P(k xmax , k ymax , ~), of an array with n stations located at (x i , y i , i=1, n) can be use to infer the location of an energy source in the wavenumber domain [for details, see Capon 1969a], as: ( 1) where ~is the angular frequency, k x and k y are wavenumbers in the N-S and E-W trends, respectively, and U ij (~) is the cross-power spectrum of i and j stations. The azimuth of the energy source to the midpoint of the array and the phase velocity of seismic waves can be calculated by arctan (k ymax /k xmax ) and /(k xmax 2 +k ymax 2 ) 1/2 , respectively.

Data analyses and results
A magnetic network with eight magnetometers was established in 1988 in Taiwan [Yen et al. 2004]. The station sites that were selected are far from populated areas and from all visible iron objects and power lines, to avoid unwanted or man-made 'noise' and interference. The geomagnetic field was continually recorded every 5 or 10 min. In 2002, three new stations were added to provide better coverage of northern to southern central Taiwan. The sampling rate of the 11 magnetometers was set for every 1 min [Chen et al. 2009b]. To further record anomalies associated with earthquakes, the sampling rates of these magnetometers were reset to 1 s in 2007. Due to one station (Kimen) located over 210 km away from the others, 10 of these stations on Taiwan Island are used in this study (Table 1). Since some of these 10 stations have occasional data gaps, the midpoint of the FK array is subject to change (see Figure  1). Figure 2 shows the longitude and latitude of the midpoint during the entire study period. It can be seen that the midpoint was relatively stationary during November 2007 to June 2008, but fluctuated a lot before and after these times. Here, the largest earthquake (EQ14, M L = 5.57; Table 2) occurred during the stationary period, and it is taken as an example to see whether seismo-magnetic signals can be recorded at the stations. We computed the amplitude at the frequency of 0.05 Hz over the 30 days before and 30 days after the earthquake occurrence (EQ14), using the Morlet wavelet transform [Farge 1992, Torrence andCompo 1998]. Note that the frequency of approximately 0.05 Hz is generally considered as the most promising frequency to detect seismo-magnetic anomalies [Fraser-Smith et al. 1990]. Due to the 1 Hz sampling rate, we derived an anomaly bond of two standard deviations, 2v, using the 86,400 amplitudes computed daily. The numbers of anomalous amplitudes were counted, as either larger or smaller than 2v. If the numbers simultaneously increased at different stations during the earthquake, seismo-magnetic signals were considered to have been recorded.
To cross-compare the stations, we standardized the numbers at each station by dividing the associated mean within 30 days before and after an earthquake. Figure 3a, b shows the variations of the standardized numbers at nine of the stations (due to the YH data gaps) 30 days before and after the April 23, 2008, M L = 5.57 earthquake (EQ14). The standardized numbers were seen to be larger for 14 days and 6 days before, as well as 8 days and 11 days after EQ14. Figure  3c shows that the Rkp index reaches 34 for 28 days before EQ14, while the standardized numbers simultaneously increased slightly. Since the Rkp on 14 days and 6 days before EQ14 are 22 and 13, respectively (both Rkp less than 24), the two days are in a magnetic quiet. It can be seen that the two anomalies intermittently increased in the standardized numbers before EQ14, which agrees with observations of discontinuous seismo-magnetic anomalous signals appearing prior to the Chi-Chi earthquake [Yen et al. 2004].
To further confirm SMEWs associated with the earthquake directions and to estimate their onset times, the azimuth and phase velocity of the geomagnetic data at the  frequency of 0.05 Hz were sequentially computed every 30 min. These computed azimuths in a moving window over 30 days with a step of 1 day are used to construct the monitored distributions. A background distribution was also developed by superimposing all of the computed azimuths during the entire period. Figure 4a shows a major peak of the PROPAGATION OF SEISMO-MAGNETIC SIGNALS  background distribution at approximately 110˚, inferring the magnetic waves at the frequency of 0.05 Hz were intrinsically westward propagations. The phase velocity of the magnetic waves was mainly around 5 ×10 4 km/s, which is slower than the velocity of light (Figure 4b). We further normalized both the background and the monitored distributions, and computed the difference proportion, which is the difference between the two distributions divided by the background distribution. A positive (or negative) value of the difference proportion means that magnetic waves propagate frequently (or rarely) at a certain azimuth. Figure 5 shows the difference proportions and the azimuths of the 26 epicenters versus the midpoint during the entire observation period. Both of the azimuths between -20˚ and 20˚ towards and away from the epicenter 10 days before the earthquakes were examined. As the azimuth and temporal resolution are 10˚ and 1 day, there are therefore 55 difference proportion values within each towards or away examined region. The averages of the difference proportion values of the two examined regions were computed, and the larger ones are indicated ( Figure 5).
To confirm that SMEWs were detected, various thresholds were tested. The marked (overall) ratios were calculated, which are the difference proportion values larger than a tested threshold divided by 55 (19800 = 36 × 550). If a marked ratio is greater than the overall ratio, we then declare that SMEWs have been detected. Figure 6 illustrates the success detection rate of the earthquakes versus the threshold. The suitable threshold was found to be 30%, which detected 62% (= 16/26) of the earthquakes. We further examined the maximum thresholds, which can successfully detect SMEWs in these 16 events, to relate these to the earthquake magnitude, but no obvious relationships could be found.

Discussion and conclusions
Seismo-magnetic anomalies have been considered as a function of the distance from the epicenter to the observation station (i.e. the epicentral distance). However, epicenters are vertically projected points from the hypocenters to the Earth surface, and they provide very limited physical information. Chen et al. [2009a] showed that seismo-magnetic fields are outwards at an interior position near the SMART, where there is an intersection between the Earth surface, the extended fault, and the auxiliary fault plane retrieved from fault-plane solutions of earthquakes [see also , Chen et al. 2010], and inwards underground at an exterior area before earthquakes. This suggests that the northern pole is close to the Earth surface and the southern pole will be underground. Geomagnetic waves propagate away from the SMART if they are outward of the northern pole and inward of the southern one. On the other hand, geomagnetic waves emitted from any northern pole can propagate to the very southern pole. As a result, geomagnetic waves propagate toward the SMART outside the seismo-magnetic fields. Similar to the aforementioned study (i.e. epicenters as reference points), we conducted the angle arrivals of SMEWs based on the SMART. Figure 5 shows that the SMART azimuths are rather more suitable than the epicenter azimuths to be related to the positive difference proportions. Although the suitable threshold is also 30%, 77% (= 20/26) of the earthquakes can be detected by using the SMART to substitute the epicenters ( Figure 6). Meanwhile, the maximum thresholds are roughly proportional to the earthquake magnitude. This means that the difference proportions that result from SMEWs are proportional to the earthquake magnitude. To further examine SMEWs away from and towards the epicenters as well as the SMART, the odds test was used, which is the quantity p/(1-p), where p is the probability of success [Agresti 2002, see also, Chen et al. 2009a]. With the odds >1, a success is more likely than a failure. Inversely, a failure is more likely than a success. Figure 7 illustrates the propagation azimuths of SMEWs versus the distances from the midpoints to the epicenters and the SMART. The odds approach 1, and therefore there is no conspicuous relationship between the propagation azimuths of the SMEWs and the distances from the midpoint to the epicenters (Figure 7a). In contrast, the propagation azimuths of the SMEWs are away from and towards the SMART where the distances from the midpoint to the SMART are <100 km and >200 km (Figure 7b), respectively. The odds Figure 6. The success detection rate versus the thresholds. Solid circles and diamonds indicate the relationships between the success detection rate and distinct thresholds using the epicenter or the SMART references, respectively. test in this study and the observation of a large decrease in the geomagnetic fields during the Chi-Chi earthquake [Yen et al. 2004], are in good agreement.
In conclusion, the azimuth and phase velocity for SMEWs carried by a scalar of the geomagnetic total intensity fields can be detected several days before an earthquake occurrence using FK analysis. On the basis that FK analysis is a standard technology for studying wave propagations, the results accordingly explain the existence of SMEWs before earthquakes. When the SMART substitutes for the epicenters to reconstruct the relationship between SMEWs and the distance to the midpoint, the success detection rate increases by 15%. The difference proportion values can be used as an indicator of the earthquake magnitude due to these positive relationships. SMEWs propagates away from and towards the SMART inside and outside the area, respectively, where they are dominated by seismogeomagnetic fields during seismogenic processes. This suggests that the SMART is a more suitable reference point than the epicenters in seismo-magnetic studies.