Time variation of the electromagnetic transfer function of the earth estimated by using wavelet transform.

In order to accurately estimate the geomagnetic transfer functions in the area of the volcano Mt. Iwate (IWT), we applied the interstation transfer function (ISTF) method to the three-component geomagnetic field data observed at Mt. Iwate station (IWT), using the Kakioka Magnetic Observatory, JMA (KAK) as remote reference station. Instead of the conventional Fourier transform, in which temporary transient noises badly degrade the accuracy of long term properties, continuous wavelet transform has been used. The accuracy of the results was as high as that of robust estimations of transfer functions obtained by the Fourier transform method. This would provide us with possibilities for routinely monitoring the transfer functions, without sophisticated statistical procedures, to detect changes in the underground electrical conductivity structure.

Introduction

It has long been suggested that the characteristics of the electrical conductivity structure of the earth may change due to the evolution of crustal activities such as seismicity, and volcanism. Attempt to verify such changes by monitoring the electric/magnetic field variations have been made since 1960s.[1)–4)] So far, however, the results of these attempts have not been really convincing due mainly to the insufficient accuracy of the obtained transfer functions (see below) and their time resolution. Traditionally transfer functions have been estimated in the frequency domain using the Fourier transform.[3),4)] The geomagnetic field data are contaminated by noises which include many transient outliers of various origins. In the conventional Fourier transform analysis, transient outliers badly affect the overall spectral properties and seriously degrade the accuracy of the transfer functions.

Since the 1980s, so called robust algorithms have been introduced mainly in the field of electromagnetic exploration.[5),6)] Some of them have demonstrated high performance for transfer function estimation. The present work is an effort to explore alternative approach to attain high accuracy of the transfer functions, through which we may be able to run routine operations to monitor their time changes.

The three components of geomagnetic field variations measured at a point can generally be expressed as; where (ω) is the frequency, and X(ω), Y(ω), and Z(ω) are the Fourier coefficients of northward, eastward, and downward components of the geomagnetic field variations and A(ω) and B(ω) which transfer the horizontal components to the vertical are called the geomagnetic transfer functions.[1)] Short-term variations in the geomagnetic field are composed of the primary source field, secondary induced field and noises of various origins. The transfer functions depend on the subterranean electrical conductivity structure which can be spatially heterogeneous and time dependent.

In dealing with the transfer functions, in this study, we applied the so called interstation transfer function (ISTF) method.[3),7),8)] In this method, transfer functions between the geomagnetic field data at a site in the area of interest and the data taken at the same time at another observation point, called reference station, are sought. Therefore, the ISTFs contain the information on the difference of the electrical structure under the two observation points, and their time variation provides the information on the relative structural change between the two points.

In this study, the station of our interest was Iwate (IWT) in the area of volcano Mt. Iwate (Fig. 1). Geomagnetic and geoelectric temporary station IWT was set up by the International Riken Frontier Research Project on Earthquakes in June 1998 as Mt. Iwate started to show signs, mainly volcanic earthquakes, of re-activation in the beginning of 1998.[9)–11)] The station was in operation until February 2002. A flux-gate magnetometer (Chiba Electronics Ltd.) was used for measuring the three components of the geomagnetic field with sampling interval of 1 s.[9)–11)]

Fig. 1.

Map showing the observation sites.[11)] IWT is the observation site for the monitoring of crustal activities. KAK indicates the Kakioka Magnetic Observatory, JMA, of which data have been used as the remote reference. Asterisk denotes the epicenter of Iwate-ken Nairiku Hokubu Earthquake (M6.2, 03 September 1998, depth 10 km).

Mt. Iwate has long been under intensive surveillance by Tohoku University and their study showed, mainly by geodetic and seismic means, that the 1998–99 activity was caused by magma intrusion to a shallow depth and its migration to the west without surface eruption.[12),13)] The activity calmed down by the late 2000. However, during this period, Iwate-ken Nairiku Hokubu Earthquake (M6.2, depth 10 km) occurred about 10 km southwest of the mountain (see also Fig. 1) on September 3, 1998.[12)] This event was unusually large for a volcanic earthquake, and its relationship to the volcanic activity is under debate.[13),14)]

The existence of pre-seismic electric/magnetic signal emission from seismic source has been postulated for some time.[15),16)] In the case of Iwate-ken Nairiku Hokubu Earthquake, it has been reported that, beyond a clear co-seismic signal,[10)] a prominent Seismic Electric Signal (SES) was recorded on August 20, i. e., two weeks before the earthquake.[17)] Moreover, the polarization (ratio between spectral intensities of vertical and horizontal components) of ULF geomagnetic variation also showed a significant increase at about the same time of SES appearance,[9),11),18)] supporting the emission from the focal zone. To find out if it was possible to see whether or not any change of transfer function caused by structural change also occurred in association with these activities was the interest of the present study. For our purpose, it was desirable that the reference station was located far enough in an area seismically quiet and sufficiently low in artificial noise during the observation period. We selected the Kakioka Magnetic Observatory, Japan Meteorological Agency as the remote reference point (Fig. 1).

To perform the ISTF, we used the wavelet transform which is known to be better fitted than Fourier transform for time-frequency analysis of noisy data due to its robustness to transient outliers.[19)] In this study we mainly follow the wavelet algorithm developed by Harada and his colleagues.[20)]

Method

ISTF:

In the geomagnetic mid-latitude regions, including Japan, the ULF geomagnetic field fluctuations originated in the upper atmosphere can be approximated by vertically incident uniform plane waves.[7)] Under this condition, the relationship of the geomagnetic field data observed at a site, and at the remote reference station may be expressed as follows[7),8),20)]; or

Here H(ω) and H(ω) are the Fourier coefficients of observed geomagnetic field data at a site and at the remote reference station. T(ω) is the interstation transfer function (ISTF), and δH(ω) denotes the uncorrelated noise. The diagonal elements, T and T, indicate the complex-valued spectral ratio in each component between the two stations. If the incident wave and the underground structure are spatially uniform, the real components of the diagonal elements are 1.0. In actual case, the diagonal elements are not exactly 1.0, due mainly to structural heterogeneity, and the non-diagonal elements, T and T, have non-zero values. The elements, T and T, have the meaning analogous to the single station transfer functions, A(ω) and B(ω) in Eq. [1].

The ISTF is usually estimated using data taken when the ratio of H(ω) to δH(ω) is high, i. e., the nighttime data on days of high solar-terrestrial geomagnetic activity. The elements of the ISTF are determined by making use of the usual power-spectrum approach.[21)]

Wavelet transform analysis:

The mother wavelet is chosen according to the signal properties and the purpose of analysis. We chose the Morlet wavelet, because it has a simple relationship between frequency and scale.[22)] The Morlet wavelet uses the exponential e as the basis function in combination with a Gaussian window function; where the dimensionless parameters t and ω0 express the time and frequency, respectively. This wavelet function meets the admissibility condition at ω0 ≥ 5.[22)] The waveform of the Morlet wavelet is illustrated in Fig. 2.

Fig. 2.

Waveform of the Morlet wavelet (ω0=5 in Eq. [4]).

Data analysis and results

We first applied the wavelet method to the IWT and KAK data for 27 hour interval from 21:00 (LT) on September 26, 1999. Fig. 3 indicates the real and imaginary parts of the resulted ISTFs as function of frequency (period). In this figure, open circles and corresponding error bars indicate the ISTFs and the upper and lower boundaries of 95% of confidence limits. The ISTFs appear quite accurate since the geomagnetic activity was extremely high on that day (Kp ≥ 5). For comparison, analysis has also been made by the algorithm which is based on the Fourier transform and robust estimation.[6)] Its results shown by solid circles and broken lines indicate remarkable coincidence with our wavelet results.

Fig. 3.

Frequency dependence of the interstation transfer function (ISTF) calculated between IWT and KAK. ISTFs were calculated for the 27 hour interval from 21:00 (LT) on September 26, 1999. Open circles and corresponding error bars indicate the obtained ISTFs and the upper and lower boundaries of 95% of confidence limits, which were calculated by using the wavelet-based approach. Solid circles and broken lines indicate the ISTFs and error bars calculated by using the sophisticated algorithm which is based on the Fourier transform.[6)]

Procedures for assessing the time variation of ISTF were as follows. We first chose the days of high solar-terrestrial geomagnetic activity (Kp ≥ 4) for every month. The number of such days was usually 3–5 per month. There were several months that had no significant storm. Night-time 6 hours (from 21:00 to 03:00 LT) of the chosen days were segmented into consecutive windows with length which would contain enough number of data points to reduce the effects of isolated outliers. The length, therefore, was set depending on the interested period of oscillation: For shorter periods, 10 to 20 cycle length was required, whereas for periods longer than, say, 1,000 s, 5 cycle length was enough.

The correlation function[23)] and the interstation transfer functions were calculated for every chosen time-frequency window. Then, the windows with interstation correlation function higher than 0.90[20)] were selected because the transfer functions obtained from those windows were expected to be higher in accuracy. The final transfer functions of the month were the results of the least-square estimate from these windows.[20)]

Fig. 4 illustrates the ISTF for the period from July 1998 to February 2002. Here only the real part of diagonal elements T and T in the periods of (a) 60 s, (b) 100 s, and (c) 1,000 s are shown. Because the intensity and occurrence frequency of geomagnetic activity were variable, the quality of the estimations also inevitably varied. Nevertheless, some of the changes seem to exceed the error bars. For example, in the year 2001, when seismic and volcanic activity was low, vague quasi-annual looking variation might be noticed in 100 s period (Fig. 4 (b)).

Fig. 4.

Time dependence of the ISTF at periods of (a) 60 s, (b) 100 s, and (c) 1,000 s for the real components of T and T. Error bars denote the upper and lower boundaries of the 95% of confidence limits. Shaded area means the period when seismic/volcanic activity of Mt. Iwate was active. (d) illustrates the daily frequency and magnitude of earthquakes at Mt. Iwate since January 1998 to March 2002 (after Monthly Report on Volcano Activities, Sendai District Meteorological Observatory, JMA). (e) illustrates the monthly integrated energy of earthquakes (log10 E=1.5M+11.8) in 30 km × 20 km area around Mt. Iwate.

Discussion

The daily frequency of earthquakes around Mt. Iwate and seismic energy release in a wider area surrounding Mt. Iwate are illustrated in Fig. 4 (d) and (e). Even though the accuracy of ISTF has been considerably improved by the present approach, no definite correlation between seismicity/volcanism and T and T can be drawn from Fig. 4. Some lowering in T and T of 60 s period in the later months of 1998 and its possible correlation with the seismic activity may be barely inferred. The present results, however, are obviously inadequate to draw any conclusion for several reasons including: a) Because of insufficient S/N ratio of available data, stable solutions were not obtained for the transfer functions in the period range shorter than 60 s. To resolve the possible structural changes due to seismic/volcanic activities in shallower depth, transfer functions for shorter periods would have been needed. b) The data before 1998–99 activity was unavailable so that the behaviors of the ISTF in normal time are unknown. c) Magnetic storms were too rare in the month of November 1998, so that no stable ISTF was obtained even for 60 s. d) As mentioned in the introduction, in the case of the Iwate-ken Nairiku Hokubu Earthquake, unusual enhancement of the geomagnetic polarization in the ULF band occurred two weeks before the earthquake. Since we used the same IWT geomagnetic data as theirs, there is a possibility that our results might have been affected by the EM emission.

Concluding remarks

So far, sensible discussion on the time changes of the transfer functions for a long term record has been difficult to make, because the conventional approaches encountered difficulties in handling noise. In this study, we have applied the wavelet transform to the calculation of transfer functions. The accuracy of the transfer functions has been improved to the level which would permit routine acquisition of their time variation by monitoring. However, due to several practical inadequacies, the present results are not convincing enough draw any definite conclusions on the long standing very important problem of existence or nonexistence of changes of subterranean electric structure related with earthquakes/magma activities. We intend to further pursue the solution of the problem by strengthening observation system, especially by adopting higher rate of data sampling, and accumulating reliable data at multiple stations.