Layout 6 Azimuthal propagation of seismo-magnetic signals from large earthquakes in Taiwan Chieh-Hung Chen1, Jann-Yenq Liu2,3,*, Tao-Ming Chang4, Ta-Kang Yeh5, Chung-Ho Wang1, Strong Wen6, Horng-Yuan Yen7, Katsumi Hattori8, Ching-Ren Lin1, Yi-Ru Chen1 1 Academia Sinica, Institute of Earth Sciences, Taipei, Taiwan 2 National Central University, Institute of Space Science, Chung-Li, Taiwan 3 National Central University, Center for Space and Remote Sensing Research, Chung, Taiwan 4 National Center for Research on Earthquake Engineering, Taipei, Taiwan 5 National Taipei University, Department of Real Estate and Built Environment, New Taipei City, Taiwan 6 National Chung Cheng University, Institute of Seismology, Chia-yi, Taiwan 7 National Central University, Institute of Geophysics, Chung-Li, Taiwan 8 Chiba University, Graduate School of Science, Inage, Chiba, Japan ANNALS OF GEOPHYSICS, 55, 1, 2012; doi: 10.4401/ag-5326 ABSTRACT 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 seismo- magnetic 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 seismo- magnetic emission waves towards and away from SMART reveal the associated anomalous propagation. 1. 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, Hayakawa et al. 1999, 2000, Gotoh et al. 2002, Hattori 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, 2007, Hattori et al. 2006]. Chen et al. [2009a] examined such magnetic anomalies statistically before 181 earthquakes (Mw ≥4.3 or ML ≥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 Special Issue: EARTHQUAKE PRECURSORS 63 Article history Received July 24, 2011; accepted November 29, 2011. Subject classification: Earthquake emission, Frequency wavenumber method, Earthquake forecast. 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 (Mw >4.3; Table 2) is also carried out, which were retrieved from the broadband array in Taiwan [Kao et al. 2002, Liang et al. 2003, 2004] for the seismology for the period from July 2007 to December 2008. 2. Theory and the FK method It has been known for many years that electromagnetic waves propagate at the speed of light of about 3 ×105 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 ×106 km (as 3 ×105/ 0.05). 64 CHEN ET AL. Figure 1. Map of Taiwan showing the locations of the 10 stations (red triangles) and 26 epicenters (blue circles). Green squares, the midpoint. Station Code Lat. (˚N) Long. (˚E) Liyutan LY 24.3467 120.7675 Tsengwen TW 23.2514 120.5167 Hengchun HC 21.935 120.8008 Yeheng YH 24.671 121.3671 Shuanlung SL 23.7902 120.9441 Pingtung PT 22.7035 120.6496 Neicheng NC 24.7181 121.6681 Hualien HL 24.0678 121.6006 Yuli YL 23.3506 121.2856 Taitung TT 22.8019 121.0519 Table 1. Geomagnetic stations in Taiwan. 65 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(kxmax, kymax, ~), of an array with n stations located at (xi, yi, 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, kx and ky are wavenumbers in the N-S and E-W trends, respectively, and Uij(~) 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 (kymax/kxmax) and �/(kxmax 2+kymax 2)1/2, respectively. 3. 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). =y x( , , ) ( ) ( ) ( )expP k k ik x x ik y yx ij i j y i j j n i n 11 1 $ - + -~ ~U == - 6 @' 1// PROPAGATION OF SEISMO-MAGNETIC SIGNALS Table 2. Details of the 26 earthquakes. Code Year Month Day Hour Min Second Lat. (˚N) Long. (˚E) Mw 1 2007 7 16 23 42 52.6 23.57 121.55 4.38 2 2007 7 23 13 40 2.1 23.72 121.64 5.07 3 2007 8 9 0 55 47.8 22.65 121.08 5.09 4 2007 8 29 1 15 29.1 21.90 121.36 4.66 5 2007 8 29 3 0 16.1 21.95 121.32 5.12 6 2007 9 6 17 51 25.2 24.28 122.25 6.17 7 2007 10 11 3 5 2.5 24.75 121.85 4.36 8 2007 10 17 14 39 59.5 23.47 121.71 4.70 9 2007 12 5 1 41 42.9 23.07 121.19 4.61 10 2008 2 17 20 33 3 23.31 121.46 5.03 11 2008 2 29 16 58 7.6 24.00 122.53 4.77 12 2008 3 4 17 31 47.8 23.21 120.70 4.89 13 2008 4 14 15 39 45.7 22.83 121.33 4.71 14 2008 4 23 18 28 42 22.87 121.68 5.57 15 2008 4 23 22 4 15.4 22.83 121.69 4.66 16 2008 5 10 19 42 1.1 23.95 122.53 5.43 17 2008 5 13 18 27 55.6 22.77 121.04 4.89 18 2008 6 1 16 59 23.9 24.86 121.79 5.12 19 2008 6 11 4 33 36.4 21.88 121.00 4.31 20 2008 8 1 18 55 49.5 24.06 121.55 4.72 21 2008 9 9 7 43 13.5 24.57 122.71 5.23 22 2008 9 10 11 55 34.1 25.08 122.22 4.37 23 2008 12 2 3 16 53.3 23.28 121.60 4.87 24 2008 12 7 21 18 36.7 23.84 122.17 4.63 25 2008 12 23 0 4 45.6 22.95 120.57 4.88 26 2008 12 30 1 31 29.1 24.66 122.36 4.39 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, ML = 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 and Compo 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, ML = 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 CHEN ET AL. 66 Figure 2. Variations in the midpoint computed by the stations from July 2007 to December 2008, according to longitude (a) (˚E) and latitude (b) (˚N). Figure 3. The standardized numbers and the Rkp index 30 days before and 30 days after EQ14. (a) Standardized numbers at NC, YH, LY and HL. (b) Standardized numbers at HC, PT, TW, TT and SL. (c) Variations un the Rkp index. 67 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 Figure 4. Distributions of (a) the derived azimuth and (b) the phase velocity. Figure 5. The difference proportion distribution. Blue to red indicate the distint difference proportions. White and yellow rectangles show the larger realtionships between the difference proportions and the azimuths of the 26 epicenters, and the SMART versus the midpoint, respectivley. The widths of the rectangles indicate a temporal period over –10 to 0 days to the earthquakes. 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 ×104 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. 4. 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 CHEN ET AL. 68 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. 69 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 seismo- geomagnetic fields during seismogenic processes. This suggests that the SMART is a more suitable reference point than the epicenters in seismo-magnetic studies. Acknowledgements. The authors wish to thank Dr. B.S. Huang and Dr. W.T. Liang of the Institute of Earth Sciences at Academia Sinica in Taiwan for providing the earthquake catalog with the fault-plane solutions from the Broadband Array in Taiwan for Seismology, as well as the Central Weather Bureau for the geomagnetic data. This study is partially supported by project NSC98-2116-M-008-006-MY3 and NSC 100-2116-M- 001-027-MY3 of the National Science Council. References Agresti, A. (2002). Categorical Data Analysis, 2 ed., John Wiley, New York, p. 710. Bernardi, A., A.C. Fraser-Smith, P.R. McGill and O.G. Vil- lard Jr. (1991). ULF magnetic field measurements near the epicenter of the Ms 7.1 Loma Prieta earthquake, Phys. Earth Planet. In., 68, 45-63. Capon, J., R.J. Greenfield and R.J. Kolter (1967). Multidi- mensional maximum-likelihood processing of a large aperture seismic array, P. IEEE, 55,192-211. Capon, J. (1969a). High-resolution frequency-wavenumber spectrum analysis, P. IEEE, 57, 1408-1419. Capon, J. (1969b). Investigation of long-period noise at the large aperture seismic array, J. Geophys. Res., 74, 3182- 3194. Chen, C.H., J.Y. Liu, W.H. Yang, H.Y. Yen, K. Hattori, C.R. Lin and Y.H. Yeh (2009a). SMART analysis of geoma- gnetic data observed in Taiwan, Phys. Chem. Earth, 34, 350-359; doi: 10.1016/j.pce.2008.09.002. Chen, C.H., C.R. Lin, H.L. Chao, H.Y. Yen, J.Y. Liu and Y.H. Yeh (2009b). Evaluation of the applicability of Chapman-Miller method on variation of the geoma- gnetic total intensity field in Taiwan from 1988 to 2007, Terr. Atmos. Ocean. Sci., 20, 799-806; doi: 10.3319/ TAO.2009.02.03.01(T). PROPAGATION OF SEISMO-MAGNETIC SIGNALS Figure 7.Relationships between the propagation azimuth of SMEWs and the distances from the midpoints to the epicenters (a) and the SMART (b) as well as their odds results. Solid circles plotted on 1 of the upper and lower vertical axes indicate that SMEWs propagate away from and towards the reference points, respectively. The shallow lines are the odds obtained with the window of five continuous earthquakes, for determination of whether the SMEW propagation azimuths are significant or insignificant over the background. Note that the maximum of the odds is 5, due to the window size. Chen, C.H., J.Y. Liu, P.Y. Lin, W.T. Liang, H.Y. Yen, K. Hat- tori and X. Zeng (2010). Pre-seismic geomagnetic ano- maly and earthquake location, Tectonophysics, 489, 240-247; doi: 10.1016/j.tecto.2010.04.018. Chen, C.H., S. Wen, J.Y. Liu, T.K. Yeh, C.H. Wang, H.Y. Yen, K. Hattori and C.R. Lin (2011). Seismomagnetic si- gnal comparison using the Morlet wavelet method, Di- saster Adv., 4, 53-60. Cheng, D.K. (1989). Field and Wave Electromagnetic, 2 ed., Addison-Wesley Pub. Comp, p. 703. Farge, M. (1992). Wavelet transforms and their applications to turbulence, Annu. Rev. Fluid Mech., 24, 395-458. Fraser-Smith, A.C., A. Bernardi, P.R. McGill, M.E. Ladd, R.A. Helliwell and O.G. Villard Jr. (1990). Low-frequency ma- gnetic field measurements near the epicenter of the Ms 7.1 Loma Prieta earthquake, Geophys. Res. Lett., 17, 1465-1468; doi: 10.1029/GL017i009p01465. Goldstein, P. and R.J. Archuleta (1991a). Deterministic fre- quency-wavenumber methods and direct measure- ments of rupture propagation during earthquakes using a dense array: Theory and methods, J. Geophys. Res., 96, 6173-615. Goldstein, P. and R.J. Archuleta (1991b). Deterministic fre- quency-wavenumber methods and direct measure- ments of rupture propagation during earthquakes using a dense array: Data analysis, J. Geophys. Res., 96, 6187- 6198. Gotoh, K., Y. Akinaga, M. Hayakawa and K. Hattori (2002). Principal component analysis of ULF geomagnetic data for Izu islands earthquakes in July 2000, J. Atmos. Elec., 22, 1-12. Hattori, K., Y. Akinaga, M. Hayakawa, K. Yumoto, T. Nagao and S. Uyeda (2002). ULF magnetic anomaly pre- ceding the 1997 Kagoshima earthquakes, In: Seismo Electromagnetics: Lithosphere-Atmosphere-Ionosphere coupling, edited by M. Hayakawa and O.A. Molchanov, Terra Scientific. Pub. Comp., Tokyo, p. 477. Hattori, K. (2004). ULF geomagnetic changes associated with large earthquakes, Terr. Atmos. Ocean. Sci., 15, 329-360. Hattori, K., I. Takahashi, C. Yoshino, N. Isezaki, H. Iwa- saki, M. Harada, K. Kawabata, E. Kopytenko, Y. Kopy- tenko, P. Maltsev, V. Korepanov, O. Molchanov, M. Hayakawa, Y. Noda, T. Nagao and S. Uyeda (2004a). ULF geomagnetic field measurements in Japan and some recent results associated with Iwateken Nairiku Hokubu earthquake in 1998, Phys. Chem. Earth, 29, 481-494. Hattori, K., A. Serita, K. Gotoh, C. Yoshino, M. Harada, N. Isezaki and M. Hayakawa (2004b). ULF geomagnetic anomaly associated with 2000 Izu islands earthquake swarm, Japan, Phys. Chem. Earth, 29, 425-435. Hattori, K., A. Serita, C. Yoshino, M. Hayakawa and N. Ise- zaki (2006). Singular spectral analysis and principal component analysis for signal discrimination of ULF geomagnetic data associated with 2000 Izu Island ear- thquake swarm, Phys. Chem. Earth, 31, 281-291. Hayakawa, M. and Y. Fujinawa (1994). Electromagnetic Phenomena Related to Earthquake Prediction, Terra Scientific. Pub. Comp., Tokyo, p. 677. Hayakawa, M. (Ed.) (1999). Atmospheric and Ionospheric Electromagnetic Phenomena Associated with Ear- thquakes, Terra Scientific. Pub. Comp., Tokyo, p. 996. Hayakawa, M., T. Itoh and N. Smirnova (1999). Fractal ana- lysis of ULF geomagnetic data associated with the Guam on August 8, 1993, Geophys. Res. Lett., 26, 2797- 2800. Hayakawa, M., T. Itoh, K. Hattori and K. Yumoto (2000). ULF electromagnetic precursors for an earthquake at Biak, Indonesia on February 17, 1996, Geophys. Res. Lett., 27, 1531-1534. Hayakawa, M. and O.A. Molchanov (2002). Seismo-Electro- magnetics: Lithosphere-Atmosphere-Ionosphere cou- pling, Terra Scientific. Pub. Comp., Tokyo, p. 477. Kao, H., Y.H. Liu, W.T. Liang and W.P. Chen (2002). Source parameters of regional earthquake in Taiwan: 1999- 2000 including the Chi-Chi earthquake sequence, Terr. Atmos. Ocean. Sci., 13, 279-298. Kawate, R., O.A. Molchanov and M. Hayakawa (1998). Ultra-low-frequency magnetic fields during the Guam earthquake of 8 August 1993 and their interpretation, Phys. Earth Planet. In., 105, 229-238. Kopytenko, Y.A., T.G. Matishvili, P.M. Voronov, E.A. Ko- pytenko and O.A. Molchanov (1993). Detection of ultra-low-frequency emissions connected with the Spi- tak earthquake and its aftershock activity, based on geo- magnetic pulsations data at Dusheti and Vardzia observatories, Phys. Earth Planet. In., 77, 85-95. Liang, W.T., Y.H. Liu and H. Kao (2003). Source parame- ters of regional earthquake in Taiwan: January–De- cember 2001, Terr. Atmos. Ocean. Sci., 14, 249-260. Liang, W.T., Y.H. Liu and H. Kao (2004). Source parame- ters of regional earthquake in Taiwan: January–De- cember 2002, Terr. Atmos. Ocean. Sci., 15, 727-741. Merzer, M. and S.L. Klemperer (1997). Modeling low-fre- quency magnetic-field precursor to the Loma Prieta earthquake with a precursory increase in fulat-zone conductivity, Pure Appl. Geophys.,150, 217-248. Molchanov, O.A., Y.A. Kopytenko, R.M. Voronov, E.A. Ko- pytenko, T.G. Matiashvili, A.C. Fraser-Smith and A. Bernardi (1992). Results of ULF magnetic field measu- rements near the epicenters of the Spitak (Ms = 6.9) and Loma Prieta (Ms = 7.1) earthquakes: Comparative analysis, Geophys. Res. Lett., 19, 1495-1498. Molchanov, O.A., M. Hayakawa and V.A. Rafalsky (1995). Penetration characteristics of electromagnetic emis- CHEN ET AL. 70 71 sions from an underground seismic source into the at- mosphere, ionosphere, and magnetosphere, J. Geophys. Res., 100, 1691-1712. Satoh, T., H. Kawase, T. Iwata, S. Higashi, T. Sato, K. Iri- kura and H.C. Huang (2001). S-wave velocity structure of the Taichung Basin, Taiwan, estimated from array and single-station records of microtremors, B. Seismol. Soc. Am., 91, 1267-1282. Schisselé, E., J. Guilbert, S. Gaffet and Y. Cansi1 (2004). Ac- curate time-frequency-wavenumber analysis to study coda waves, Geophys. J. Int., 158, 577-591; doi: 10.1111/ j.1365-246X.2004.02211.x. Telesca, L., M. Balasco, G. Colangelo, V. Lapenna and M. Macchiato (2004). Investigating the multifractal pro- perties of geoelectrical signals measured in southern Italy, Phys. Chem. Earth, 29, 295-303. Telesca, L. and K. Hattori (2007). Nonuniform scaling be- havior in ultra low frequency (ULF) earthquake-related geomagnetic signals, Physica A, 384, 522-528; doi: 10.1016/j.phys a.2007.05.040. Torrence, C. and G.P. Compo (1998). A practical guide to wavelet analysis, B. Am. Meteorol. Soc., 79, 61-78. Yen, H.Y., C.H. Chen, Y.H. Yeh, J.Y. Liu, C.R. Lin and Y.B. Tasi (2004). Geomagnetic fluctuations during the 1999 Chi-Chi earthquake in Taiwan, Earth Planets Space, 56, 39-45. *Corresponding author: Jann-Yenq Liu, National Central University, Institute of Space Science, Chung-Li, Taiwan email: jyliu@jupiter.ss.ncu.edu.tw. © 2012 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. PROPAGATION OF SEISMO-MAGNETIC SIGNALS