- Email: [email protected]
ULF geomagnetic changes possibly associated with the 2008 Iwate–Miyagi Nairiku earthquake
Journal of Asian Earth Sciences 41 (2011) 442–449
Contents lists available at ScienceDirect
Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jseaes
ULF geomagnetic changes possibly associated with the 2008 Iwate–Miyagi Nairiku earthquake Takuya Hirano, Katsumi Hattori * Graduate School of Science, Chiba University, 1-33, Yayoi, Inage, Chiba 263-8522, Japan
a r t i c l e
i n f o
Article history: Available online 1 June 2010 Keywords: ULF geomagnetic anomaly Spectrum density ratio analysis Wavelet transform The 2008 Iwate–Miyagi Nairiku earthquake Normalized spectrum density ratio
a b s t r a c t There are many reports on earthquake-related electromagnetic phenomena. Anomalous ULF geomagnetic field changes associated with earthquake is one of the most convincing and promising phenomena due to deeper skin depth. Since ULF signals associated with large earthquakes are weak, effective signal discrimination methods should be required. Several methods for the signal discrimination have been developed so far: which are spectrum density ratio analysis, geomagnetic transfer function analysis, fractal analysis, principal component analysis, direction finding analysis, and so on. In this study, we investigate ULF geomagnetic changes possibly associated with the 2008 Iwate–Miyagi Nairiku earthquake based on spectral density ratio analysis. Geomagnetic data observed at Esashi, where the epicentral distance is about 47 km and Kakioka, the distance is about 317 km, and as a reference station have been analyzed. Wavelet transform have been performed for the spectral density analysis instead of the conventional FFT method. Before the earthquake, the variation of spectral density ratio, Sz/Sx and Sz/Sy, at the nearest station of Esashi exhibits an apparent increase from the trend. On the contrary, there are no corresponding significant changes at a remote station of Kakioka. After investigating the singularity of the increase using normalized spectrum density ratio, the enhancement is the most significant in intensity and duration for the all analyzed period. The level of peak is beyond the 3r and its duration is 3 days. The lead time is about 3– 4 weeks before the earthquake. At the periods from 3 to 105 s, similar anomalous changes occurred. These facts suggest the anomalous change is a possible candidate of earthquake-related ULF magnetic change. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Recently, there are many reports on seismo-electromagnetic phenomena in a wide frequency range (Hayakawa and Molchanov, 2002; Molchanov and Hayakawa, 2008). Measurements of electromagnetic phenomena can be classified into three types: (1) passive groundbased observation, (2) active ground-based observation using transmitter signals, and (3) satellite observation (e.g. Hattori, 2004). Among of these observational methods, one of the most promising methods is earthquake-related ULF (ultra low frequency, with frequency less than 1 Hz) phenomena because of skin depth. ULF geomagnetic data observed at the ground are superposition of several signals. The most intense signals are geomagnetic pulsations. The next intense signals are artificial noises due to DC-driven train, factories and so on. The third signals are earthquake-related signals. Earthquake related signals are very weak, so there have been developed several kind of methods for
* Corresponding author. Tel.: +81 43 290 2801; fax: +81 43 290 2859. E-mail address: [email protected] (K. Hattori). 1367-9120/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jseaes.2010.04.038
discriminating these signals from other signals: (1) spectrum density ratio analysis by means of the ratio of vertical component to horizontal component (Hayakawa et al. 1996; Hattori et al., 2002, 2004; Saroso et al., 2009), (2) geomagnetic transfer function analysis (Yanagihara and Nagano, 1976; Saroso et al., 2009), (3) monoand multi-fractal analyses (Hayakwa et al., 1999; Ida et al., 2005; Telesca et al., 2008), (4) principal component analysis (Hattori et al., 2004, 2006), (5) direction finding analysis (Ismaguilov et al., 2003), and so on. Fig. 1 represents the relationship between magnitude and epicentral distance for earthquake-related ULF geomagnetic anomalies (Hattori, 2004). The detectable distance for anomaly is given by the formula 0.025R < M 4.5, where R is epicentral distance and M is magnitude. In this paper, ULF geomagnetic data observed at Esashi station are analyzed, where the epicentral distance is about 47 km. Therefore, it is expected to detect some ULF anomalies before the day of the 2008 Iwate–Miyagi Nairiku earthquake. In this paper, the ULF geomagnetic changes associated with the 2008 Iwate–Miyagi Nairiku earthquake are presented based on the spectrum density ratio analysis using wavelet transform.
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
443
Fig. 1. The relationship between magnitude and epicentral distance for ULF anomaly (after Hattori, 2004).
2. The 2008 Iwate–Miyagi Nairiku earthquake and ULF geomagnetic stations The 2008 Iwate–Miyagi Nairiku earthquake occurred at N39.03°, E140.89° on June 13, 2008 with magnitude M7.2 and the depth is 8 km. The source mechanism is reverse faulting with (Strike, Dip, Rake) = (203°, 37°, 93°) (Japan Meteorological Agency, 2008; Earthquake Research Institute, 2008). The earthquake triggered landslide. Over 10 people died and about 450 people injured. Furthermore, several micro-earthquakes (M < 1.4) occurred about 40 min before the earthquake near the epicenter (National Research Institute for Earth Science and Disaster Prevention, 2008) and three aftershocks with M > 5 occurred on June 14, 2008 (M5.7 and M5.2) and June 16, 2008 (M5.3) (Japan Meteorological Agency, 2008). Fig. 2 shows a map of the epicenter and ULF geomagnetic stations. The data observed at Esashi and Kakioka are used. The distances from the epicenter are 47 km and 317 km, respectively. At the nearest Esashi station, induction type 3 component magnetometer is in operation by Geographical Survey Institute, Japan (Satoh et al., 2004). The sampling frequency is 15 Hz. In this study, the data are re-sampled to 1 Hz. As a reference, Kakioka data are adopted, where a three component fluxgate type magnetometer is in operation with 1 Hz sampling by Japan Meteorological Agency. We analyze 3 years data from January 1, 2006 to 3, December 2008. Fig. 2. The map of epicenter and observation sites.
3. Data analysis 3.1. Wavelet spectrum analysis In this paper, spectral density analysis with wavelet transform is performed for detection of anomalous changes in observed data instead of conventional FFT. Because FFT analysis is not suitable for transient signals such as impulses. The concept of this spectral density analysis is based on the idea that signals associated with
solar-terrestrial origins are likely to be plane waves with vertically incident. Therefore, signals originated from solar-terrestrial have no vertical component of magnetic field theoretically. If the Esashi observatory is placed over an uniform 1D Earth structure or over 2D electrical conductivity anomaly (on its central part), vertical component of induction field becomes zero. However, the vertical component has some values because of heterogeneity of the underground structure. Therefore, the ratio of horizontal compo-
444
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
Fig. 3. Wavelet sonograms at Esashi station on April 1, 2008. (a) North–south component, (b) east–west component, and (c) vertical component.
nent to vertical component seems to become a certain value. On the other hand, crustal origin signal from the additional underground current is considered to have a lager vertical component except for just above the current position. This suggests that the spectral density ratio such as Sz/Sx and Sz/Sy shows the higher value in the active period of earthquake-related ULF electromagnetic phenomena. Here Sx, Sy, and Sz mean the spectral density at a certain frequency for north–south, east–west, and vertical components, respectively. In this paper, Morlet wavelet (Morlet et al., 1982) is selected for a mother function. Morlet wavelet is defined by w0 ðgÞ with a sinusoidal function and a Gaussian window as shown
w0 ðgÞ ¼ p1=4 eibg eg
2 =2
ð1Þ
where g is a non-dimensional time parameter. b is the non-dimensional frequency and here taken to be 6 to satisfy the admissibility condition (Farge, 1992). The continuous wavelet transform of a time series data xn is defined as the convolution of xn with a scaled and translated version of w0 ðgÞ as shown
0 ðn nÞdt W n ðsÞ ¼ xn w s n0 ¼0 N1 X
ð2Þ
where * indicates the complex conjugate. s and n are the wavelet scale and the localized time index n, respectively. The relation between wavelet scale s and Fourier frequency f is given as follows:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 0 b þ 2 þ b2 1 A S¼@ f 4p
ð3Þ
As for the reconstruction of the time series data, we use the sum of the real part of the wavelet transform over all scales (Farge, 1992),
xn ¼
j djdt 1=2 X Re W n ðSj Þ 1=2 C d w0 ð0Þ j¼0 Sj
ð4Þ 1=2
The factor w0 ð0Þ removes the energy scaling, while the Sj converts the wavelet transform to an energy density. w0 ð0Þ and C d are constant of w1=4 and 0.776 for our Morlet wavelet, respectively (Torrence and Compo, 1998). 3.2. Spectrum density analysis Figs. 3 and 4 show examples of wavelet sonograms observed at Esashi station on April 1, 2008 (normal day) and on May 19, 2008 (abnormal day (25 days before the main shock)), respectively. It is found the intensity of the vertical component below 100 s increase in Fig. 4c. According to the previous reports, the possible anomalous changes are dominant in range of 10–100 s. The results seem to be consistent. In order to investigate tendency of variation in
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
445
Fig. 4. Wavelet sonograms at Esashi station on May 19, 2008. (a) North–south component, (b) east–west component, and (c) vertical component.
power, 3 h averaging is performed. The results in the period of 10.41 s in 2008 are given in Fig. 5a–c corresponds to the variation of Sx, Sy, and Sz in dB, respectively. Fig. 5d indicates the variation of RKp, which represents geomagnetic activity caused by solar-terrestrial interaction. When RKp > 30, the geomagnetic activity is considered active. Fig. 5e shows the variation of Es index, which is the daily sum of local earthquake energy in terms of magnitude in the following formula (Hattori et al., 2006).
Es ¼ log
X 104:8þ1:5M r2 1day
ð5Þ
where M and r indicate magnitude of earthquake and hypocenter distance from Esashi station. In this paper, earthquakes with r < 200 km are selected from the earthquake catalog issued by Japan Meteorological Agency. The day of the mainshock is indicated by a vertical red line. In Fig. 5, among many spikes, the pulse immediately before the earthquake shows the maximum intensity in 2008 and the duration seems to be longer than other spikes. At least there is no significant geomagnetic storm during the period of the anomalous enhancement appearance. However, we cannot con-
sider, with confidence, that it is a distinguished anomaly. Similar tendency is observed in other periods. In order to investigate subterranean effects more clearly, the variation of spectral density ratio such as Sz/Sx and Sz/Sy is taken. Fig. 6 indicates the variation of averaged spectral density ratios over 1 day, Sz/Sx and Sz/Sy, at Esashi station for the period of 10.41 s together with RKp and Es index in comparison with Fig. 5.1 The thin blue and thick red lines in each panel indicate the 1 day value and its ±5 days running mean, respectively. It is found that spectrum density ratios increase apparently 35–20 days before the earthquake. There is no significant geomagnetic disturbance around the day of the earthquake. Fig. 7 shows the variation of averaged spectral density ratios over 1 day, Sz/Sy, at: (a) Esashi station and (b) Kakioka reference station for the period of 10.41 s together with RKp and Es index for a whole analyzed period. There is a data missing for Esashi station from day 251 to 270. The variation of spectrum density ratio at
1 For interpretation of color in Figs.3–8, the reader is referred to the web version of this article.
446
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
Fig. 5. The variation of power spectrum at the period of 10.41 s at Esashi station in 2008. (a) North–south component, (b) east–west component, (c) vertical component, and (d) RKp index, and (e) Es index.
Fig. 6. The variation of spectrum density ratio at the period of 10.41 s in 2008. (a) Sz/Sx, (b) Sz/Sy, (c) RKp index, and (d) Es index.
Esashi reveals an enhancement before the earthquake from the baseline as shown in Fig. 6. However, there is no anomalous enhancement for the ratio at Kakioka before the earthquake. The similar tendency is obtained in the Sz/Sx (and also for the periods of 3–33 s). These facts are highly suggestive of a local or an anomalous change at Esashi. However, there are some increases of the ratio through the whole analyzed period. In order to eliminate the seasonal variation and examine singularity, the following normalization is performed:
SZ SZ =Si SZ =Si ¼ ; Si r
i ¼ x; y
ð6Þ
where Sz/Si is ±5 days running average, Sz =Si is ±30 days running average and r is the standard deviation for 30 days data. We call Sz Si
normalized spectrum density ratio hereafter and it is defined
the threshold for anomaly using r. In this paper, the criterion of anomaly is set at +2.5r. Fig. 8 illustrates variations of normalized spectrum density ratio Sz/Sy at periods of: (a) 10.41 s, (b) 52.48 s, and (c) 104.95 s for a whole period at Esashi. From Fig. 8a, there exist two anomalies in the period of 10.41 s and the anomaly 29–27 days before the earthquake has the highest value beyond 3 r and the longest duration (3 days). Similar characteristics are obtained in the period of 52.48 s and 101.95 s. In fact, the anomaly appear in the period from
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
447
Fig. 7. The variation of spectrum density ratio (Sz/Sy) at the period of 10.41 s. (a) Esashi station, (b) Kakioka station, (c) RKp index, (d) Es index at Esashi, and (e) Es index at Kakioka.
Fig. 8. The variations of normalized spectrum density ratio at Esashi. (a) Period of 10.41 s, (b) period of 52.48 s, (c) period of 104.95 s, (d) RKp index, and (e) Es index.
448
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
Fig. 9. The relationship between magnitude and epicentral distance for ULF anomaly. The results of this study are given by a white star (Esashi) and a dark star (Kakioka) (after Hattori (2004)).
3 s to 105 s, although the range of analyzed periods is from 3 s to 1060 s. The anomaly before the earthquake is a wide-band frequency phenomenon. These facts are highly suggestive of a local or an anomalous change at Esashi and a possible candidate of earthquake-related ULF change. Possible mechanisms to generate the underground current system are: (1) microfracturing (Molchanov and Hayakawa, 1998) and (2) electrokinetic effect (Ishido and Muzutani, 1981). Theoretically, electromagnetic signals derived from these models have wideband spectrum. The earth crust acts as a low pass filter. Therefore, it is supposed that longer period signals can be detected on the ground.
In summary, it is important to use multiple station data for the spectral density ratio analysis of ULF magnetic data for monitoring crustal activity. Especially to remove the global effect, the usage of an adequate reference site is essential. Moreover, the analysis of the longer data is also important for estimate the normal variation at the site. The proposed methodology of the normalized spectral density ratio analysis provides one of the statistical approaches for this aim. However, the number of the ULF geomagnetic station is small at the moment to investigate further study on the spatial distribution of the anomalies. Furthermore, in order to make convince the result, other methodologies such as geomagnetic transfer function analysis and fractal analysis should be performed. Acknowledgements
4. Conclusion The ULF geomagnetic data observed at Esashi and Kakioka have been analyzed to detect any possible precursory signatures associated with the 2008 Iwate–Miyagi Nairiku earthquake. Before the earthquake, the variation of spectral density ratio, Sz/Sx and Sz/Sy, at the nearest station of Esashi exhibits an apparent increase from the trend. On the contrary, there are no corresponding significant changes at a remote station of Kakioka. After investigating the singularity of the increase using normalized spectrum density ratio, the enhancement is the most significant in intensity and duration for the all analyzed period. The level of peak is beyond the 3r and its duration is 3 days. The lead time is about 3–4 weeks before the earthquake. At the periods from 3 to 105 s, similar anomalous changes occurred. These facts suggest the anomalous change is a possible candidate of earthquake-related ULF magnetic change. Therefore the Fig. 1 can be updated as Fig. 9. Our study has added the two data points, and the empirical threshold seems to be still valid.
The authors would like to express thanks to Geographical Survey Institute and Kakioka Magnetic Observatory, Japan Meteorological Agency for providing geomagnetic data at Esashi and Kakioka. Also we thank to Japan Meteorological Agency for the earthquake data. This research was partly supported by a Grandin-Aid for Scientific Research of Japan Society for Promotion of Science (19403002) and National Institute of Information and Communication Technology (R&D promotion funding international joint research). Reference Earthquake Research Institute, 2008. The University of Tokyo Iwate–Miyagi Nairiku Earthquake Preliminary Finite Fault Model.. Farge, M., 1992. Wavelet transforms and their applications to turbulence. Annual Review of Fluid Mechanics 24 (1), 395–457. Hattori, K., 2004. ULF geomagnetic changes associated with large earthquakes. Terrestrial, Atmospheric and Oceanic Sciences 15, 329–360.
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449 Hattori, K., Akinaga, Y., Hayakawa, M., Yumoto, K., Nagao, T., Uyeda, S., 2002. ULF magnetic anomaly preceding the 1997 Kagoshima earthquakes. In: Hayakawa, M., Molchanov, O.A. (Eds.), Seismo Electromagnetics: Lithosphere– Atmosphere–Ionosphere Coupling. TERRAPUB, Tokyo, pp. 19–28. Hattori, K., Takahashi, I., Yoshino, C., Isezaki, N., Iwasaki, H., Harada, M., Kawabata, K., Kopytenko, E., Kopytenko, Y., Maltsev, P., Korepanov, V., Molchanov, O., Hayakawa, M., Noda, Y., Nagao, T., Uyeda, S., 2004. ULF geomagnetic field measurements in Japan and some recent results associated with Iwateken Nairiku Hokubu earthquake in 1998. Physics and Chemistry of the Earth 29 (4–9), 481–494. Hattori, K., Serita, A., Yoshino, C., Hayakawa, M., Isezaki, N., 2006. Singular spectral analysis and principal component analysis for signal discrimination of ULF geomagnetic data associated with 2000 Izu Island Earthquake Swarm. Physics and Chemistry of the Earth 31 (4–9), 281–291. Hayakawa, M., Ito, T., Smirnova, N., 1999. Fractal analysis of ULF geomagnetic data associated with the Guam earthquake on August 8, 1993. Geophysical Research Letters 26 (18), 2797–2800. Hayakawa, M., Kawate, R., Molchanov, O.A., Yumoto, K., 1996. Results of ultralowfrequency magnetic field measurements during the Guam earthquake of 8 August, 1993. Geophysical Research Letters 23, 241–244. Hayakawa, M., Molchanov, O.A. (Eds.), 2002. Seismo Electromagnetics Lithosphere– Atmosphere–Ionosphere Coupling. TERAPUB, Tokyo. Ida, Y., Hayakawa, M., Adalev, A., Gotoh, K., 2005. Multifractal analysis for the ULF geomagnetic data during the 1993 Guam earthquake. Nonlinear Processes in Geophysics 12 (2), 157–162. Ismaguilov, V.S., Kopytenko, Yu.A., Hattori, K., Hayakawa, M., 2003. Variations of phase velocity and gradient values of ULF geomagnetic disturbances connected with the Izu strong earthquakes. Natural Hazards and Earth System Science 3 (3–4), 211–215. Japan Meteorological Agency, 2008. (in Japanese).
449
Ishido, T., Muzutani, J., 1981. Experimental and theoretical basis of electrokinetic phenomena in rock–water systems and its applications to geophysics. Journal of Geophysical Research 86 (B3), 1763–1775. Molchanov, O.A., Hayakawa, M., 1998. On the generation mechanism of ULF seismogenic electromagnetic emissions. Physics of the Earth and Planetary Interiors 105 (3–4), 201–210. Molchanov, O.A., Hayakawa, M., 2008. Seismo-electromagnetics and Related Phenomena: History and Results. TERRAPUB, Tokyo. 189p. Morlet, J., Arens, G., Fourgeau, E., Giard, D., 1982. Wave propagation and sampling theory – part I: complex signal and scattering in multilayered media. Geophysics 47, 203–221. National Research Institute for Earth Science and Disaster Prevention, 2008.
Contents lists available at ScienceDirect
Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jseaes
ULF geomagnetic changes possibly associated with the 2008 Iwate–Miyagi Nairiku earthquake Takuya Hirano, Katsumi Hattori * Graduate School of Science, Chiba University, 1-33, Yayoi, Inage, Chiba 263-8522, Japan
a r t i c l e
i n f o
Article history: Available online 1 June 2010 Keywords: ULF geomagnetic anomaly Spectrum density ratio analysis Wavelet transform The 2008 Iwate–Miyagi Nairiku earthquake Normalized spectrum density ratio
a b s t r a c t There are many reports on earthquake-related electromagnetic phenomena. Anomalous ULF geomagnetic field changes associated with earthquake is one of the most convincing and promising phenomena due to deeper skin depth. Since ULF signals associated with large earthquakes are weak, effective signal discrimination methods should be required. Several methods for the signal discrimination have been developed so far: which are spectrum density ratio analysis, geomagnetic transfer function analysis, fractal analysis, principal component analysis, direction finding analysis, and so on. In this study, we investigate ULF geomagnetic changes possibly associated with the 2008 Iwate–Miyagi Nairiku earthquake based on spectral density ratio analysis. Geomagnetic data observed at Esashi, where the epicentral distance is about 47 km and Kakioka, the distance is about 317 km, and as a reference station have been analyzed. Wavelet transform have been performed for the spectral density analysis instead of the conventional FFT method. Before the earthquake, the variation of spectral density ratio, Sz/Sx and Sz/Sy, at the nearest station of Esashi exhibits an apparent increase from the trend. On the contrary, there are no corresponding significant changes at a remote station of Kakioka. After investigating the singularity of the increase using normalized spectrum density ratio, the enhancement is the most significant in intensity and duration for the all analyzed period. The level of peak is beyond the 3r and its duration is 3 days. The lead time is about 3– 4 weeks before the earthquake. At the periods from 3 to 105 s, similar anomalous changes occurred. These facts suggest the anomalous change is a possible candidate of earthquake-related ULF magnetic change. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Recently, there are many reports on seismo-electromagnetic phenomena in a wide frequency range (Hayakawa and Molchanov, 2002; Molchanov and Hayakawa, 2008). Measurements of electromagnetic phenomena can be classified into three types: (1) passive groundbased observation, (2) active ground-based observation using transmitter signals, and (3) satellite observation (e.g. Hattori, 2004). Among of these observational methods, one of the most promising methods is earthquake-related ULF (ultra low frequency, with frequency less than 1 Hz) phenomena because of skin depth. ULF geomagnetic data observed at the ground are superposition of several signals. The most intense signals are geomagnetic pulsations. The next intense signals are artificial noises due to DC-driven train, factories and so on. The third signals are earthquake-related signals. Earthquake related signals are very weak, so there have been developed several kind of methods for
* Corresponding author. Tel.: +81 43 290 2801; fax: +81 43 290 2859. E-mail address: [email protected] (K. Hattori). 1367-9120/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jseaes.2010.04.038
discriminating these signals from other signals: (1) spectrum density ratio analysis by means of the ratio of vertical component to horizontal component (Hayakawa et al. 1996; Hattori et al., 2002, 2004; Saroso et al., 2009), (2) geomagnetic transfer function analysis (Yanagihara and Nagano, 1976; Saroso et al., 2009), (3) monoand multi-fractal analyses (Hayakwa et al., 1999; Ida et al., 2005; Telesca et al., 2008), (4) principal component analysis (Hattori et al., 2004, 2006), (5) direction finding analysis (Ismaguilov et al., 2003), and so on. Fig. 1 represents the relationship between magnitude and epicentral distance for earthquake-related ULF geomagnetic anomalies (Hattori, 2004). The detectable distance for anomaly is given by the formula 0.025R < M 4.5, where R is epicentral distance and M is magnitude. In this paper, ULF geomagnetic data observed at Esashi station are analyzed, where the epicentral distance is about 47 km. Therefore, it is expected to detect some ULF anomalies before the day of the 2008 Iwate–Miyagi Nairiku earthquake. In this paper, the ULF geomagnetic changes associated with the 2008 Iwate–Miyagi Nairiku earthquake are presented based on the spectrum density ratio analysis using wavelet transform.
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
443
Fig. 1. The relationship between magnitude and epicentral distance for ULF anomaly (after Hattori, 2004).
2. The 2008 Iwate–Miyagi Nairiku earthquake and ULF geomagnetic stations The 2008 Iwate–Miyagi Nairiku earthquake occurred at N39.03°, E140.89° on June 13, 2008 with magnitude M7.2 and the depth is 8 km. The source mechanism is reverse faulting with (Strike, Dip, Rake) = (203°, 37°, 93°) (Japan Meteorological Agency, 2008; Earthquake Research Institute, 2008). The earthquake triggered landslide. Over 10 people died and about 450 people injured. Furthermore, several micro-earthquakes (M < 1.4) occurred about 40 min before the earthquake near the epicenter (National Research Institute for Earth Science and Disaster Prevention, 2008) and three aftershocks with M > 5 occurred on June 14, 2008 (M5.7 and M5.2) and June 16, 2008 (M5.3) (Japan Meteorological Agency, 2008). Fig. 2 shows a map of the epicenter and ULF geomagnetic stations. The data observed at Esashi and Kakioka are used. The distances from the epicenter are 47 km and 317 km, respectively. At the nearest Esashi station, induction type 3 component magnetometer is in operation by Geographical Survey Institute, Japan (Satoh et al., 2004). The sampling frequency is 15 Hz. In this study, the data are re-sampled to 1 Hz. As a reference, Kakioka data are adopted, where a three component fluxgate type magnetometer is in operation with 1 Hz sampling by Japan Meteorological Agency. We analyze 3 years data from January 1, 2006 to 3, December 2008. Fig. 2. The map of epicenter and observation sites.
3. Data analysis 3.1. Wavelet spectrum analysis In this paper, spectral density analysis with wavelet transform is performed for detection of anomalous changes in observed data instead of conventional FFT. Because FFT analysis is not suitable for transient signals such as impulses. The concept of this spectral density analysis is based on the idea that signals associated with
solar-terrestrial origins are likely to be plane waves with vertically incident. Therefore, signals originated from solar-terrestrial have no vertical component of magnetic field theoretically. If the Esashi observatory is placed over an uniform 1D Earth structure or over 2D electrical conductivity anomaly (on its central part), vertical component of induction field becomes zero. However, the vertical component has some values because of heterogeneity of the underground structure. Therefore, the ratio of horizontal compo-
444
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
Fig. 3. Wavelet sonograms at Esashi station on April 1, 2008. (a) North–south component, (b) east–west component, and (c) vertical component.
nent to vertical component seems to become a certain value. On the other hand, crustal origin signal from the additional underground current is considered to have a lager vertical component except for just above the current position. This suggests that the spectral density ratio such as Sz/Sx and Sz/Sy shows the higher value in the active period of earthquake-related ULF electromagnetic phenomena. Here Sx, Sy, and Sz mean the spectral density at a certain frequency for north–south, east–west, and vertical components, respectively. In this paper, Morlet wavelet (Morlet et al., 1982) is selected for a mother function. Morlet wavelet is defined by w0 ðgÞ with a sinusoidal function and a Gaussian window as shown
w0 ðgÞ ¼ p1=4 eibg eg
2 =2
ð1Þ
where g is a non-dimensional time parameter. b is the non-dimensional frequency and here taken to be 6 to satisfy the admissibility condition (Farge, 1992). The continuous wavelet transform of a time series data xn is defined as the convolution of xn with a scaled and translated version of w0 ðgÞ as shown
0 ðn nÞdt W n ðsÞ ¼ xn w s n0 ¼0 N1 X
ð2Þ
where * indicates the complex conjugate. s and n are the wavelet scale and the localized time index n, respectively. The relation between wavelet scale s and Fourier frequency f is given as follows:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 0 b þ 2 þ b2 1 A S¼@ f 4p
ð3Þ
As for the reconstruction of the time series data, we use the sum of the real part of the wavelet transform over all scales (Farge, 1992),
xn ¼
j djdt 1=2 X Re W n ðSj Þ 1=2 C d w0 ð0Þ j¼0 Sj
ð4Þ 1=2
The factor w0 ð0Þ removes the energy scaling, while the Sj converts the wavelet transform to an energy density. w0 ð0Þ and C d are constant of w1=4 and 0.776 for our Morlet wavelet, respectively (Torrence and Compo, 1998). 3.2. Spectrum density analysis Figs. 3 and 4 show examples of wavelet sonograms observed at Esashi station on April 1, 2008 (normal day) and on May 19, 2008 (abnormal day (25 days before the main shock)), respectively. It is found the intensity of the vertical component below 100 s increase in Fig. 4c. According to the previous reports, the possible anomalous changes are dominant in range of 10–100 s. The results seem to be consistent. In order to investigate tendency of variation in
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
445
Fig. 4. Wavelet sonograms at Esashi station on May 19, 2008. (a) North–south component, (b) east–west component, and (c) vertical component.
power, 3 h averaging is performed. The results in the period of 10.41 s in 2008 are given in Fig. 5a–c corresponds to the variation of Sx, Sy, and Sz in dB, respectively. Fig. 5d indicates the variation of RKp, which represents geomagnetic activity caused by solar-terrestrial interaction. When RKp > 30, the geomagnetic activity is considered active. Fig. 5e shows the variation of Es index, which is the daily sum of local earthquake energy in terms of magnitude in the following formula (Hattori et al., 2006).
Es ¼ log
X 104:8þ1:5M r2 1day
ð5Þ
where M and r indicate magnitude of earthquake and hypocenter distance from Esashi station. In this paper, earthquakes with r < 200 km are selected from the earthquake catalog issued by Japan Meteorological Agency. The day of the mainshock is indicated by a vertical red line. In Fig. 5, among many spikes, the pulse immediately before the earthquake shows the maximum intensity in 2008 and the duration seems to be longer than other spikes. At least there is no significant geomagnetic storm during the period of the anomalous enhancement appearance. However, we cannot con-
sider, with confidence, that it is a distinguished anomaly. Similar tendency is observed in other periods. In order to investigate subterranean effects more clearly, the variation of spectral density ratio such as Sz/Sx and Sz/Sy is taken. Fig. 6 indicates the variation of averaged spectral density ratios over 1 day, Sz/Sx and Sz/Sy, at Esashi station for the period of 10.41 s together with RKp and Es index in comparison with Fig. 5.1 The thin blue and thick red lines in each panel indicate the 1 day value and its ±5 days running mean, respectively. It is found that spectrum density ratios increase apparently 35–20 days before the earthquake. There is no significant geomagnetic disturbance around the day of the earthquake. Fig. 7 shows the variation of averaged spectral density ratios over 1 day, Sz/Sy, at: (a) Esashi station and (b) Kakioka reference station for the period of 10.41 s together with RKp and Es index for a whole analyzed period. There is a data missing for Esashi station from day 251 to 270. The variation of spectrum density ratio at
1 For interpretation of color in Figs.3–8, the reader is referred to the web version of this article.
446
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
Fig. 5. The variation of power spectrum at the period of 10.41 s at Esashi station in 2008. (a) North–south component, (b) east–west component, (c) vertical component, and (d) RKp index, and (e) Es index.
Fig. 6. The variation of spectrum density ratio at the period of 10.41 s in 2008. (a) Sz/Sx, (b) Sz/Sy, (c) RKp index, and (d) Es index.
Esashi reveals an enhancement before the earthquake from the baseline as shown in Fig. 6. However, there is no anomalous enhancement for the ratio at Kakioka before the earthquake. The similar tendency is obtained in the Sz/Sx (and also for the periods of 3–33 s). These facts are highly suggestive of a local or an anomalous change at Esashi. However, there are some increases of the ratio through the whole analyzed period. In order to eliminate the seasonal variation and examine singularity, the following normalization is performed:
SZ SZ =Si SZ =Si ¼ ; Si r
i ¼ x; y
ð6Þ
where Sz/Si is ±5 days running average, Sz =Si is ±30 days running average and r is the standard deviation for 30 days data. We call Sz Si
normalized spectrum density ratio hereafter and it is defined
the threshold for anomaly using r. In this paper, the criterion of anomaly is set at +2.5r. Fig. 8 illustrates variations of normalized spectrum density ratio Sz/Sy at periods of: (a) 10.41 s, (b) 52.48 s, and (c) 104.95 s for a whole period at Esashi. From Fig. 8a, there exist two anomalies in the period of 10.41 s and the anomaly 29–27 days before the earthquake has the highest value beyond 3 r and the longest duration (3 days). Similar characteristics are obtained in the period of 52.48 s and 101.95 s. In fact, the anomaly appear in the period from
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
447
Fig. 7. The variation of spectrum density ratio (Sz/Sy) at the period of 10.41 s. (a) Esashi station, (b) Kakioka station, (c) RKp index, (d) Es index at Esashi, and (e) Es index at Kakioka.
Fig. 8. The variations of normalized spectrum density ratio at Esashi. (a) Period of 10.41 s, (b) period of 52.48 s, (c) period of 104.95 s, (d) RKp index, and (e) Es index.
448
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449
Fig. 9. The relationship between magnitude and epicentral distance for ULF anomaly. The results of this study are given by a white star (Esashi) and a dark star (Kakioka) (after Hattori (2004)).
3 s to 105 s, although the range of analyzed periods is from 3 s to 1060 s. The anomaly before the earthquake is a wide-band frequency phenomenon. These facts are highly suggestive of a local or an anomalous change at Esashi and a possible candidate of earthquake-related ULF change. Possible mechanisms to generate the underground current system are: (1) microfracturing (Molchanov and Hayakawa, 1998) and (2) electrokinetic effect (Ishido and Muzutani, 1981). Theoretically, electromagnetic signals derived from these models have wideband spectrum. The earth crust acts as a low pass filter. Therefore, it is supposed that longer period signals can be detected on the ground.
In summary, it is important to use multiple station data for the spectral density ratio analysis of ULF magnetic data for monitoring crustal activity. Especially to remove the global effect, the usage of an adequate reference site is essential. Moreover, the analysis of the longer data is also important for estimate the normal variation at the site. The proposed methodology of the normalized spectral density ratio analysis provides one of the statistical approaches for this aim. However, the number of the ULF geomagnetic station is small at the moment to investigate further study on the spatial distribution of the anomalies. Furthermore, in order to make convince the result, other methodologies such as geomagnetic transfer function analysis and fractal analysis should be performed. Acknowledgements
4. Conclusion The ULF geomagnetic data observed at Esashi and Kakioka have been analyzed to detect any possible precursory signatures associated with the 2008 Iwate–Miyagi Nairiku earthquake. Before the earthquake, the variation of spectral density ratio, Sz/Sx and Sz/Sy, at the nearest station of Esashi exhibits an apparent increase from the trend. On the contrary, there are no corresponding significant changes at a remote station of Kakioka. After investigating the singularity of the increase using normalized spectrum density ratio, the enhancement is the most significant in intensity and duration for the all analyzed period. The level of peak is beyond the 3r and its duration is 3 days. The lead time is about 3–4 weeks before the earthquake. At the periods from 3 to 105 s, similar anomalous changes occurred. These facts suggest the anomalous change is a possible candidate of earthquake-related ULF magnetic change. Therefore the Fig. 1 can be updated as Fig. 9. Our study has added the two data points, and the empirical threshold seems to be still valid.
The authors would like to express thanks to Geographical Survey Institute and Kakioka Magnetic Observatory, Japan Meteorological Agency for providing geomagnetic data at Esashi and Kakioka. Also we thank to Japan Meteorological Agency for the earthquake data. This research was partly supported by a Grandin-Aid for Scientific Research of Japan Society for Promotion of Science (19403002) and National Institute of Information and Communication Technology (R&D promotion funding international joint research). Reference Earthquake Research Institute, 2008. The University of Tokyo Iwate–Miyagi Nairiku Earthquake Preliminary Finite Fault Model.
T. Hirano, K. Hattori / Journal of Asian Earth Sciences 41 (2011) 442–449 Hattori, K., Akinaga, Y., Hayakawa, M., Yumoto, K., Nagao, T., Uyeda, S., 2002. ULF magnetic anomaly preceding the 1997 Kagoshima earthquakes. In: Hayakawa, M., Molchanov, O.A. (Eds.), Seismo Electromagnetics: Lithosphere– Atmosphere–Ionosphere Coupling. TERRAPUB, Tokyo, pp. 19–28. Hattori, K., Takahashi, I., Yoshino, C., Isezaki, N., Iwasaki, H., Harada, M., Kawabata, K., Kopytenko, E., Kopytenko, Y., Maltsev, P., Korepanov, V., Molchanov, O., Hayakawa, M., Noda, Y., Nagao, T., Uyeda, S., 2004. ULF geomagnetic field measurements in Japan and some recent results associated with Iwateken Nairiku Hokubu earthquake in 1998. Physics and Chemistry of the Earth 29 (4–9), 481–494. Hattori, K., Serita, A., Yoshino, C., Hayakawa, M., Isezaki, N., 2006. Singular spectral analysis and principal component analysis for signal discrimination of ULF geomagnetic data associated with 2000 Izu Island Earthquake Swarm. Physics and Chemistry of the Earth 31 (4–9), 281–291. Hayakawa, M., Ito, T., Smirnova, N., 1999. Fractal analysis of ULF geomagnetic data associated with the Guam earthquake on August 8, 1993. Geophysical Research Letters 26 (18), 2797–2800. Hayakawa, M., Kawate, R., Molchanov, O.A., Yumoto, K., 1996. Results of ultralowfrequency magnetic field measurements during the Guam earthquake of 8 August, 1993. Geophysical Research Letters 23, 241–244. Hayakawa, M., Molchanov, O.A. (Eds.), 2002. Seismo Electromagnetics Lithosphere– Atmosphere–Ionosphere Coupling. TERAPUB, Tokyo. Ida, Y., Hayakawa, M., Adalev, A., Gotoh, K., 2005. Multifractal analysis for the ULF geomagnetic data during the 1993 Guam earthquake. Nonlinear Processes in Geophysics 12 (2), 157–162. Ismaguilov, V.S., Kopytenko, Yu.A., Hattori, K., Hayakawa, M., 2003. Variations of phase velocity and gradient values of ULF geomagnetic disturbances connected with the Izu strong earthquakes. Natural Hazards and Earth System Science 3 (3–4), 211–215. Japan Meteorological Agency, 2008.
449
Ishido, T., Muzutani, J., 1981. Experimental and theoretical basis of electrokinetic phenomena in rock–water systems and its applications to geophysics. Journal of Geophysical Research 86 (B3), 1763–1775. Molchanov, O.A., Hayakawa, M., 1998. On the generation mechanism of ULF seismogenic electromagnetic emissions. Physics of the Earth and Planetary Interiors 105 (3–4), 201–210. Molchanov, O.A., Hayakawa, M., 2008. Seismo-electromagnetics and Related Phenomena: History and Results. TERRAPUB, Tokyo. 189p. Morlet, J., Arens, G., Fourgeau, E., Giard, D., 1982. Wave propagation and sampling theory – part I: complex signal and scattering in multilayered media. Geophysics 47, 203–221. National Research Institute for Earth Science and Disaster Prevention, 2008.