Variations of Geomagnetic Response Functions before the 2011 Tohoku Earthquake ()
1. Introduction
First of all, we will explain our methodology in the analyses of geomagnetic field variations, and we will describe our research purpose on finding out any precursors to the 2011 Tohoku earthquake (EQ) from geomagnetic data with special reference to induction vectors.
1.1. Explanation of Our Analysis Methods
Geomagnetic field
(where ex, ey, and ez are unit vectors directed to North, East and downward correspondingly) continuously exists in and around the Earth and is recorded nowadays digitally with a time reading interval Δt. After the conventional processing using Fourier transform a B record of duration t2-t1 is transformed into a superposition of harmonic components with periods
.
Response function (RF) is the term widely used in natural sciences and mathematics. In the geoelectromagnetic studies of electrical conductivity σ(x, y, z) of the Earth’s interior [1] [2] [3] the RFs are supposed to be some functions derived from the Earth’s electromagnetic (EM) data that provide us with a possibility to determine the conductivity structure of the Earth. EM RFs are usually frequency/(period T) dependent and then they are complex functions having real (index u) and imaginary (index v) parts. We use two of these functions.
Induction vector
(A and B are determined from the linear equation:
). Real induction vectors
possess an important property: in the notation of Wiese used here, they are directed away from a good conductor.
Anomalous horizontal magnetic variation tensor [M] is determined from the linear system of equations
,
, where r0 and ri are coordinates of base (reference) and some other observation place. Tensor [M] reflects the change in geoelectric structures between two places.
The processing of observed geomagnetic field B for the monitoring of geodynamic and other environmental processes implies transformation of 3 component geomagnetic field time-series with time reading interval Δt (1 min or 1 s in our study) into a variety of time-series with temporal resolution Δτ (
) of different RF components: Re and Im, x and y at the set of periods
of received harmonics (
).
The theory of geoelectromagnetic methods is developed for natural source field in the form of vertically incident plane wave (Tikhonov-Cagniard (T-C) model), which usually holds for an external source field of magnetosphere-ionosphere origin (named as magnetotelluric (MT) field) for the periods less than 104 s. Ideally RF depends only on the Earth’s conductivity distribution which is sensitive to the stress variations and therefore to geodynamic processes including the earthquake (EQ) preparation.
Really the observed B(T) is composed of several sources: 1) MT-field, part of which meets the T-C model requirements, the rest does not meet them and makes an error, 2) ever-present noise, both man-made and natural and sometimes BLE—internal fields of Lithosphere Emission.
1.2. Introduction to Our Present Study
RFs and their variations, especially in relation with the occurrence of EQs and EQ preparation, were studied in Japan for many years and were described in many works between which we cite only few [4] [5] [6] [7] . However, in the last two decades this RF approach became less frequently used for EQ studies because of stronger noise at Japanese observatories due to the enhanced human activity.
It is recently agreed that different kinds of phenomena do take place prior to an EQ, including surface deformation, geochemical anomalies, EM radiation, atmospheric and ionospheric perturbations, etc. [8] . Especially, the geomagnetic and ULF (ultra low frequency) magnetic field variations among EM precursors are considered to be a promising candidate for EQ prediction, because significant precursors have been accumulated for several huge EQs [3] . Further some precursors have already been found to be statistically correlated with EQs, such as ionospheric perturbations [3] [8] .
A huge Tohoku EQ (with magnitude of 9) happened offshore of the Tohoku area on 11 March, 2011, and many scientists have tried to find any EQ precursors to this EQ. Different kinds of precursors (middle- and short-term precursors) to this 2011 Tohoku EQ have been reported so far including EM, geochemical anomalies and others, though most were recognized only retrospectively [8] . By paying particular attention again to the RF methods, we have analyzed the available geomagnetic data all over Japan in order to obtain any middle- and short-term EQ precursors in geomagnetic variations for this Tohoku EQ. Some results were already presented in Russian [9] [10] [11] , and the latest results based on the RF methods will be presented in this paper. Finally, our RF results will be compared with the corresponding former results by other researchers in the conclusion part.
2. Data Used and Data Processing
We have obtained the geomagnetic data from 17 Japanese geomagnetic observatories (see Figure 1, Table 1) with temporal resolution of 1 min. Table 1 summarizes names, codes, and geographic (and geomagnetic) coordinates of Japanese geomagnetic observatories used in this paper, and the last column indicates the period of our processing. After the processing we obtain the values of 4 components of induction vector Au, Bu, Av and Bv for each day for 5 period intervals centered at 225, 450, 1800 and 3600 s. In order to reduce a great scatter of data, every daily values were monthly averaged.
Then for the conversion of geomagnetic field B time series into RF time series we used the advanced multi-windows rr (remote reference) robust programs with coherency control [12] [13] .
![]()
Figure 1. Map of Japan with real Cu and imaginary Cv induction vectors for the period 450 s at 17 geomagnetic observatories by 1-min records data. a) Mean vectors for the year of 2001, stars present EQs with M > 7.8 since 2001, elements of plate tectonics: white arrows represent the directions of plates motions, E. P.—Eurasian Plate, P. P.—Pacific Plate, P. S. P.—Philippines Sea Plate, S. T.—Sagami Trough, dotted line—volcanic front, b) Mean vectors for 2011, in circles EQs with M ≥ 7. The depth of hypocenters is less than 50 km for all EQs.
![]()
Table 1. Geomagnetic observatories used in this paper. Station name, code, geographic (and geomagnetic) coordinates.
3. Results of Processing for the Separation of Middle-Term Precursors
Analyzing large material of processed data for 15 years from 2001 to 2015, we found that aperiodic variations (or enhancement of annual variation) of induction vectors were observed at periods of 225, 450 and 900 s during 3 - 5 years before the Tohoku EQ at the following stations of HAR, KAK, OTA, KNZ and TTK, most clearly at period 450 s presented in Figure 2. We should emphasize that no such aperiodic variations were observed at other stations including ESA and MIZ, which are the nearest to the EQ epicenter. The best correlation of middle-term anomalies is observed at the two most remote (620 km) stations HAR and TTK: at both stations we notice strong synchronous variations of induction vectors with maxima in the end of 2008, the end of 2009, with several maximums in 2010 and in the beginning of 2011, and recover to the previous level after the Tohoku EQ. We may suppose that these aperiodic variations can be the middle-term precursors of the Tohoku EQ. These observatories are located not at the shortest distance from the EQ, that is in agreement with
![]()
Figure 2. Variations of the monthly mean induction vector components at the period of 450 s during 2001-2015 at 5 observatories with significant changes before the Tohoku EQ: HAR, KAK, OTA, KNZ and TTK. Vertical bars present the uncertainty of every monthly mean value. Two strong EQs are shown by vertical lines.
well-known phenomenon of spatial selectivity of EQ precursors known during the centuries for hydrological precursors and recently proven for BLE registered in the form of Seismic Electric Signal (SES) [14] [15] [16] [17] .
Having 1 min time series we can analyze only geomagnetic variations with period T > 3 min and the most interesting shorter part of ULF spectra (0.01 - 10 Hz) where the strongest emissions (BLE) have been observed [18] [19] , is left not resolved. So, we get 1 s data for observatories KAK, KNZ, ESA, MIZ and for short time intervals for UCU (Uchiura) and KYS (Kiyosumi) in the Boso area to analyze RFs for periods T > 5 s.
4. Boso Conductivity Anomaly
Processing of records from 18 observatories (16 of them are shown in Figure 1 and KYS, UCU in Figure 3) to determine the horizontal tensors [M] with KAK as the base (reference) station yields the absence of noticeable horizontal tensor anomalies in ESA, MIZ, HAR, TTK and MUR, but reveals their existence in KNZ, UCU, OTA and KYS (Figure 3(a)). In KNZ and OTA the enhancement of real tensor components Mxx and Myy equals to ≈40% and ≈30% correspondingly at periods T < 500 s with decreasing at longer periods. This result was supported by direct visual measurements described below. At closely located observatories KNZ and OTA, the latitudinal (E-W) component of induction vector at period 450 s and shorter increased (in 2011 comparatively to 2001) in opposite directions: westward in KNZ and eastward in OTA (see Figure 1(b)). It means that between these two observatories an additional current (of geodynamic origin) appeared in 2011.
4.1. Visual Analysis of Geomagnetic Records
Considering the time series of geomagnetic field synchronous records ( [11] , Figure 2), we noticed that MT field appears synchronously at all observatories while noise appears locally at each one. The stations most contaminated by noise, are UCU and KNZ which are the nearest to DC railways. But during the after-mid-night time interval from ≈1:30 to ≈4:30 LT (16:30 - 19:30 UT) the strong noise from DC railways almost disappears, although at some other stations (HAR and partly TTK) weaker noise continues all night.
Direct measurements of the strong MT variation amplitudes in each component provide a check (not precise but very reliable) of the results obtained by our processing. So, the enhancement of Bx at KNZ and OTA at approximately 30% - 40% exists and it can be interpreted only by the electrical conductivity anomaly under the observatories, i.e. under the central part of the Boso peninsula.
4.2. Comparison with Geology and Tectonic Evidence
The relation between Mxx and Myy anomalies in KNZ defines the WNW-ESE strike of the Boso conductivity anomaly, and geological data [20] [21] presented in Figure 3(b) and Figure 3(c) suggest the existence of an anomalous conductor of WNW-ESE strike in Miura Group sediments of the Kanto plain at depth 0 - 4 km. Relations between Mxx and Myy in UCU, OTA and KYS are different as seen in Figure 3(a). It means that the direction of anomalous currents is also different under each observatory. Calculations show that at least 50% of anomalous currents should be located near the surface in the sediments of the Kanto plain to match with the received data.
![]()
Figure 3. RFs and Neogene sediments in the Kanto plain and the Boso peninsula. a) Frequency characteristics of the modulus of horizontal tensor [M] main components at Boso observatories with reference to the base observatory KAK. Every curve is a mean value for 7 years (KNZ), 25 days (UCU), 1 year (OTA) and 12 days (KYS). Realization length Δt and interval of averaging t are written at every graph (they were chosen depending on the interval of available data and their discreteness); b) Thickness of Neogene sediments [20] [21] and induction vectors for period ≈1 hour, named observatories with real and imaginary vectors—our processing, the other 6 real vectors are taken from [4] ; c) Thickness of Miura Group deposits with “N.8 half graben fills” and induction vectors obtained from the dominant DC noise and directed to the noise source—DC railways given by thin lines for suburb railways, by thick line—for magistral one; d) Similar results for other 4 periods; e) Results of the same data processing with an attempt to eliminate the effect of the noise by means of either the use of only 4 after-midnight hours or by the application of remote reference techniques; f) Frequency characteristics of induction vector components at KNZ. Dark lines for the wide periods range-ss (single station) processing of 3 days realizations with sliding reference to middle day, light lines-processing of after-midnight 4 hours records of the same 20 days from February 18 to March 9, 2011.
On the other hand, the plate tectonic evidences that the Boso anomaly is located over the Sagami trough, structure at the depth 15 - 20 km in the complex junction of three lithosphere plates (Figure 1(a)). The strike of the trough is the same, WNW-ESE, so some part of the anomalous conductor can be located in the Sagami trough.
The eastern part of both conductors (shallow sediments and deep trough) has a contact with sea water, while the western one suppositionally contacts with a magma reservoir. In such a circuit it can be some unstable area(s) with conductivity strongly dependent on stress and sensitive to stress changes related with EQs.
Figure 3(b) shows the vectors for a period of about an hour (4000 s) at which industrial noise is practically absent and those vectors adequately reflect the heterogeneity of geoelectric structure. In Figure 3(c) vectors at the period 25 s are built with dominated noise field, which is much greater than MT field at four observatories considered. Real and imaginary vectors at periods 16, 25, 50 and partly at 200 s (Figure 3(d)) are directed to the source of noises—the nearest railway. To reduce influence of the noise, nighttime records and remote reference technique were used (Figure 3(f). Received corrected vectors appeared still very scattered as in Figure 3(e) and can be used cautiously for monitoring of geodynamic processes. However, those vectors averaged over a long period of time can be used for clarifying the geoelectric structure. Real vector in KNZ at the shortest periods is found to direct to north. It means that the most conductive part of Boso electrical conductivity anomaly is located to the south of KNZ, apparently near the southern side of the asymmetric “N.8 half graben fills” of the Miura Group sediments [21] .
Results of the processing of single station and nighttime records at KNZ are given in Figure 3(f). We see that full day and nighttime results are significantly different from each other only for Bu and partly Bv components, because a railway is located to the west of KNZ and brings the distortions mainly in the eastern component. The northern component is less affected by noise at periods 100 s and more that opens the possibility to use it for the separation of EQ precursors.
In the time intervals with high MT field the RF corresponds to the conductivity distribution and can be used for the study of precursors arisen from changes of conductivity. As was discussed above the anomalous conductor under KNZ can be a part of the large scale electrical circuit (sea water—magma reservoir) strongly sensitive to stress changes. The shallow geoelectric structure of Kanto plain sediments is the part of this circuit. The changes of electrical currents in the circuit before an EQ can be recorded even with the use of strong dominant DC noise field. This idea is supported by the use of strong noise from DC railways for the study of geoelectrical structure [22] .
5. Dependence on Geomagnetic Activity
From the analysis of magnetograms and frequency characteristics of induction vector and from the study of the correlation coefficient (CC) between induction vector variations and variations of Kp index, it has been proven [10] that there is an inherent interval of periods with “quasi-parity” between noise and the MT field at each station. During the low geomagnetic activity noise is dominant, while during the high activity MT field is dominant and a RF strongly depends on geomagnetic activity. It should be emphasized that even in this case CC does not exceed 0.6 as a rule, leaving the space for separating precursors, especially using an adequate modern processing program filtering out noisy intervals.
6. Short-Term EQ Precursor Separation
6.1. Background of the Method
The induction vector derived from very noisy records, practically from noise field, has a small stable phase. Therefore, if some other magnetic field which is not so stable (let it be a precursor field), is superimposed on the field of noise, exactly the phase of induction vector will be the most sensitive component for such a precursor separation. The nature of this “other” magnetic field can be BLE (data of [23] evidence that after February 22, 2011 at an observatory UCU (in the south of the Boso peninsula) BLE appeared with maximum at periods 30 - 100 s). An alternative explanation of the precursor is the change of anomalous current in the large scale circuit including the Boso anomaly. In any case phase in such a noisy environment can be the most sensitive parameter for precursor separation.
Below we apply a new approach developed by V.I.Tregubenko [24] , who used it for processing the data of seismo-prognosis monitoring network in Ukraine. He separated precursors before few strongest (M ≈ 4) Crimean EQs occurred during 15 years, in particular before the Sudak, Crimea EQ M3.9 on 24.01.2005 [24] . We applied this approach to KNZ, KAK and ESA 1-s data, but the precursor was found only in KNZ.
6.2. Processing of KNZ Data
The processing was made with the use of the program developed by Varentsov [12] . Coherences were used as weights estimates for averaging the results. To minimize the effect of noise the estimates with multiple coherence less than 0.7 were ignored that allowed us to obtain minimally shifted estimates of induction vector’s components. The maximum anomalous effect before the Tohoku EQ was observed in the phase of the induction vector northern component—arg (A) for periods of 100 - 200 s. For the longer periods the anomaly gradually decreases, so that we now present the results for T = 100 s. The processing and analysis were performed in two steps:
Step 1. 7 years long (2005-2011) geomagnetic field time series with every 1 s reading was processed for every month as a single unit. arg(A) time series with every month reading were received as in Figure 4(a) and analyzed by a polynomial approximation approach. The most significant first, second and seventh (quasi-annual) order polynomials were extracted from the rough data and we obtained Figure 4(b) which is more perspective for comparison with EQs. But 1 month temporal resolution of RF is not sufficient for such an analysis. As for the Tohoku EQ precursors, we find a 9 month long negative anomaly in 2010 approximately 1 year before the main event.
Step 2. 2 years long (2010-2011) 1 s time series were processed for every day as single unit. A large scatter of every day results was reduced by averaging with moving window of 5 days long with 1 day shift. From these curves, i.e. arg(A) time series with every day reading, we extracted first, second and seventh (quasi-annual) order polynomials determined in Step 1. The result is shown by the grey rough curve in Figure 5(a). The bold solid line is the result of averaging with 29 days window with 1 day shift to suppress monthly variations arising from the Moon rotation around the Earth (gravity effect) [25] and the Sun rotation around its axis (magnetic activity effect). The mean of annual variation during 7 years of 2005-2011 has been subtracted, but we note a residual annual variation in Figure 5(a) (RF annual variation enhancement was recorded at several observatories during 2 - 3 years before Tohoku EQ and can be considered as middle-term precursor). So, such a residual annual variation was determined once more from 2 years data presented in Figure 5(a) and subtracted, the result is presented in Figure 5(b).
![]()
Figure 4. Variations of arg(A) for the period T = 100 s at the observatory KNZ during the period of 2005-2011. a) Monthly values of arg(A) are given by dots. Dashed straight line and bold solid line are the first- and second-order polynomials respectively. Thin smooth line is quasi-annual variations obtained as seventh-order polynomial approximation. In the upper left corner the determination of mean annual variation is shown (first and second order polynomials were excluded). b) Variations of arg(A) after removal of the first, second and seventh (quasi-annual) order polynomials. Moments of strong EQs M > 5 are indicated by vertical lines and magnitude M (in parentheses, the first number is the depth of the hypocenter, the second one—distance from epicenters to the observatory KNZ, both are given in km).
![]()
Figure 5. Variations of arg(A) for the period T = 100 s at the observatory KNZ during 2010-2011: a) After removal trends (first and second polynomial approximation) presented in Figure 4. b) The same after the removal residual annual variation determined from Figure 4(a). Rough gray curves are the results of every day processing after averaging with moving 5 days width window with 1 day shift. Bold solid lines are the result of averaging (for elimination of monthly variation) with 29 days window with 1 day shift. The occurrence time of the Tohoku EQ on 11.03.2011 is marked by a vertical arrow.
6.3. Discussion
The variations of arg(A) given by the bold solid line in Figure 5(b) are cleaned from periodic variations of annual and monthly periods, that make it more convenient to identify EQ precursors. Indeed, we see a strong bay-like variation that begins almost 2 months before the Tohoku EQ and has approximately 2 month duration. This is in good agreement with the finding of Varotsos team [26] [27] . In particular, Varotsos et al. [26] showed that the initiation of an SES activity in Japan appears almost simultaneously with the minimum of the fluctuations of the order parameter of seismicity analyzed in natural time, and such a minimum was clearly observed on January 5, 2011, i.e. almost two months before the Tohoku EQ occurrence [27] . Moreover, the spatio-temporal study of the aforementioned minimum gave a ground for a reliable estimation of the epicentral area of the impending main shock [27] .
The time of beginning of a bay-like variation and its duration depends on the magnitude of an expected EQ. This time is equal to approximately 2 weeks for the processed Crimean EQs with magnitude approximately 4 [24] . But usually the precursors are not so smooth, but they are rather some sequence of separate episodes. And we obtained such precursors a few weeks as well as some days before the Tohoku EQ by the high resolution processing [9] . Really, averaging with 29 days window eliminated not only the monthly variation but also many shorter ones. The rough grey graph (also every day values with only 5 days averaging) in Figure 5(b) can be used for a more detailed separation of precursors, but the impact of superposed background from the noise and magnetic activity should be kept in mind and eliminated as our future work.
7. Summary and Conclusion
We have calculated induction vectors using data from Japanese observatories for many years preceding the 2011 Tohoku huge EQ, and we will summarize the main results obtained in this paper as follows.
1) Thorough extensive analyses of very noisy geomagnetic data of Japanese observatories allow us to separate EQ anomalies which can be attributed to the medium-term EQ precursors to the disastrous Tohoku EQ on 11 March 2011.
2) Under the central part of Boso peninsula at a distance of 30 km from Tokyo, we could reliably detect a conductor anomaly. A change of anomalous field over the anomaly before the EQ was identified. Tectonic interpretation of the results contains a supposition that the Boso anomaly can be a fragment of large-scale anomaly between the Pacific water and the magma reservoir 50 - 100 m deep.
3) The EQ precursor can occur by perturbations of the field of strong noise from DC railway, which may be the case for the Boso anomaly. The Boso anomaly is a very attractive subject for our further detailed study, probably as a sensitive place for monitoring EQ precursors.
4) A short-term EQ precursor (two months long) just two months before the EQ was clearly identified in the phase data (less sensitive to noise) at the observatory, KNZ with the use of the method developed by Tregubenko and often used in Ukraine.
5) Finally the RF approach is found to be still a useful tool for presentation and storage of EM monitoring data even in any noisy environment.
Finally we try to compare our findings with the similar previous analyses for the same Tohoku EQ by other researchers. Kopytenko et al. [23] investigated the magnetic field variations at three stations (ESA, MIZ, and KAK as used in this paper) over the 11-year period 2000-2011. They found a medium-term anomaly beginning around 3 years before the EQ, and a short-term precursor in the frequency range of 0.033 - 0.01 Hz was observed starting from 22 February, 2011. Another type of anomalous daily variation of ULF magnetic field data at ESA was found to occur approximately 2 months prior to the EQ [28] . The timing of these anomalies as based on the B field analysis seems to be consistent with our results, but those anomalies are observed in the data of the station ESA closest to the EQ epicenter, that looks as a sharp contrast with our finding based on RF method. We think that results based on the B field analysis seem to require a further study at all observatories, taking in mind that the anomalies in geomagnetic and ULF data are usually not so intense. Hayakawa et al. [29] have performed a critical, so-called natural time analysis on the ULF data at KAK for this EQ, and found that the horizontal magnetic field (f = 0.03 - 0.05 Hz) fulfilled all criticality conditions for 3-5 March, 2011, a few days before the EQ, lending a support to our previous result with the induction vector analysis [9] .
Acknowledgements
The authors are grateful to Geospatial Information Authority of Japan, Japan Meteorological Agency and World Data Center for Geomagnetism (Kyoto) for providing us with good quality data. We are also thankful to V. I. Tregubenko for his valuable consultations on the application of his method, Drs. T. A. Klymkovych and I. M. Varentsov for their processing programs. Finally the authors would like to express their sincere thanks to Prof. K. Hattori of Chiba University for kindly providing the data in the Boso area.