ULF/ELF Atmospheric Radiation in Possible Association to the 2011 Tohoku Earthquake as Observed in China

Recordings of a search-coil magnetometer located in Yongsheng in the Yunnan province of China have been used to find atmospheric ULF/ELF precursors to the 2011 Tohoku earthquake (EQ) with M=9 which occurred on 11 March 2011 with a large epicentral distance of about 4100 km. The combined characteristics of the horizontal magnetic fields, the ratio of tangential and radial field components in the numerator and root mean square of ellipticity in the denominator, were used to detect seismo-related ULF/ELF signals. This value is found to increase before the EQ, and then to decrease during a few weeks, with a reliable maximum on 8 March (3 days before the EQ). We also analyzed azimuth distributions of pulsed radiation in the frequency range 2-5 Hz during the previous week before the main shock. This radiation reached a maximal value on 10 March one day before the EQ date, and the azimuth of its source is estimated to coincide with the position of the main shock. The results are further compared with the ULF/ELF precursors of the same EQ observed at 3 stations in Japan by Ohta et al. (2013).


Introduction
The existence of electromagnetic anomalies before large earthquakes (EQs) have been confirmed by numerous studies during the last few decades (e.g., Molchanov & Hayakawa, 2008;Hayakawa (Ed), 2012, 2013;Hayakawa, 2015;Sorokin et al., 2015).According to those previous works, it is found that there are two different kinds of electromagnetic short-term anomalies: the first category is direct effect of the lithosphere in the EQ hypocenter (or epicenter) (Molchanov & Hayakawa, 2008), and the second category is an indirect effect such as seismo-atmospheric and -ionospheric signatures (abnormal pre-EQ modification of the atmosphere and ionosphere).The first category includes DC seismic electric signals (Varotsos, 2005), ULF (ultra low frequency) radiation (Hayakawa et al., 2007a;Fraser-Smith, 2009;Kopytenko et al., 2009;Hattori, 2013), and so on.While we have to mention the publication of very few papers criticizing the presence of seismogenic ULF radiation (e.g., Campbell, 2009).As for the second category, we can list over-the-horizon VHF propagation (Hayakawa et al., 2007b;Devi et al., 2012) and seismo-atmospheric ULF/ELF (extremely low frequency) radiation (Schekotov et al., 2007;2013) as the seismo-atmospheric signature and the perturbations both in the lower ionosphere by means of subionospheric VLF (very low frequency) signals (Hayakawa et al., 2010) and as the depression of ULF magnetic field (Schekotov et al., 2006;2013) and the upper ionosphere based on the ionosonde and in-situ data as the seismo-ionospheric signatures (Liu, 2009;Parrot, 2013).Though the mechanism of those seismo-atmospheric and -ionospheric effects is not well understood, further studies are required on the mechanism of lithosphere-atmosphere-ionosphere coupling (Molchanov et al., 2004;Molchanov & Hayakawa, 2008;Freund, 2009;Pulinets & Ouzounov, 2011).
Among the above anomalies, several of them have already been found to be statistically correlated with EQs.Such examples are increasing, including the seismo-ionospheric ULF depression (Schekotov et al., 2006;2013), the seismo-atmospheric ULF/ELF radiation (Schekotov et al., 2007;2013), the ionospheric perturbations both in the lower ionosphere (by means of subionospheric VLF/LF (low frequency) transmitter signals) (Hayakawa et al., 2010) and also in the upper F region ionosphere (by means of ionosonde observations etc.) (Liu, 2009).
In addition to the above statistical studies, case studies for any huge EQs are still of great importance in understanding short-term seismo-electromagnetic effects.Again, we take the 2011 Japan EQ in this paper.Concerning this EQ occurred on 11 March, 2011, some anomalies have already been reported by different methods.The most convincing one is the lower ionospheric perturbation by Hayakawa et al. (2012a), who have examined the characteristics of VLF/LF signal transmitted from the NLK US transmitter based on the measurements at three stations in Japan, and have found a serious decrease in the night-time average amplitude of the VLF/LF signal on 5 and 6 March.A further detailed study on VLF/LF has followed (Hayakawa et al., 2013), and they have compared those ionospheric perturbations with the corresponding crustal movements (Kamiyama et al., 2014).In support to those results, Zhou et al. (2013) have found the Schumann resonance anomalies on 8 March, 3 days prior to the main shock, based on the ELF observations at two stations in China.The anomaly in Schumann resonances was interpreted in terms of the seismo-lower ionospheric perturbations (Hayakawa et al., 2005;Nickolaenko et al., 2006).At the same time, the short-time TEC (total electron content) anomaly above the EQ epicenter and its magnetically conjugate region in the ionosphere was observed to be significantly enhanced on 8 March based on the measurements of the global navigation satellite systems and ionosonde stations near the epicenter (Yao et al., 2012;Le et al., 2013).On the other hand, there are few observational evidence belonging to the first category of direct radiation from the lithosphere.Ohta et al. (2013) have examined the ULF/ELF measurements at three stations in Japan and have discovered the ULF/ELF precursor from the atmosphere as a seismo-atmospheric signature of EQs, which occurred on 6 March, about 5 days prior to the main shock.
As a further extension of our recent paper by Ohta et al. (2013), the purpose of this paper is to try to re-confirm the presence of seismogenic atmospheric ULF/ELF radiation as already observed in Japan by Ohta et al. (2013), by examining the similar ULF/ELF records by a magnetometer in China with a large epicentral distance of about 4100 km.By using the same method as in Schekotov et al. (2007) based on the examination of polarization characteristics of atmospheric ULF/ELF pulses for the same 2011 Tohoku EQ, we try to re-confirm that such seismogenic ULF/ELF radiation could be detected even in China prior to the EQ, and also we want to demonstrate the effectiveness of the analysis method developed by Schekotov et al. (2007) and used by Ohta et al. (2013).

ULF/ELF Data
An ELF station is working in a cave near the town of Yongsheng (abbreviated as YS) (geographic coordinates: 26.7°N, 100.8°E) in the Yunnan province of China.Three orthogonal search-coil magnetometers with sampling rate 100 Hz are installed, and the working frequency band is about 1-30 Hz, and the sensitivity at a frequency of 10 Hz is about 100 fT/sqrt(Hz).Figure 1 illustrates the frequency response of the amplitude of our transfer function of this magnetometer.In this work, we analyze only two horizontal magnetic field components, and use recordings during an interval covering 5 months from 1 January to 31 May 2011 including the 2011 March 11 EQ date.
Typical daily evolutions of horizontal field components in time and frequency domains on a particular quiet day of 5 March 2011 are shown in Figure 2, from which moderate level of industrial interferences can be seen.Figure 3 shows a map covering the region from China to Japan including the YS station and the epicenter of the 2011 Tohoku EQ.Herein the position of magnetometer is shown by a black diamond and marked by a letter YS.EQs with magnitude (M) >6 are plotted by colorful circles, and the seismic data are obtained from the ANSS Worldwide Earthquake Catalog (http://www.ncedc.org/anss/catalog-search.html).The color and dimension of a circle depend on the focal depth and the magnitude (M) of the EQ respectively, and the scales for depth and magnitude are shown by means of two subplots on the right bottom part of the map.The largest circle absolutely corresponds to the 2011 Tohoku EQ.

Seismo-atmospheric ULF/ELF electromagnetic radiation and signal processing methods
As is given in our first paper by Schekotov et al. (2007), seismo-atmospheric radiation in the ULF/ELF frequency band seems to provide us with a possibility of predicting an EQ; that is, not only predicting the occurrence time of a forthcoming EQ, but also predicting the direction to its source of radiation.Some later works (Schekotov et al., 2008(Schekotov et al., , 2013;;Molchanov & Hayakawa, 2008;Hayakawa et al., 2012b;Ohta et al., 2013) have confirmed that the direction of radiation source is coincident approximately with the position of the epicenter of a future EQ.The details of the analysis method can be found in Schekotov et al. (2007) and Ohta et al. (2013), but we repeat only the essential points.

Direction Finding
We select impulsive signals (with wideband spectra) whose amplitudes are above five times the average amplitude of ELF pulses, and we determine the direction of the source of such impulsive ULF/ELF radiation as being perpendicular to the main axis of polarization ellipse (See the corresponding equation in Ohta et al. (2013)).
An example showing the azimuth distribution is shown in Figure 4.The azimuthal distribution of emissions is laid over the map of region of our interest shown in equidistant azimuthal projection.The distribution of azimuthal angle (α) is represented by an angle histogram, which is a polar plot showing the distribution of α values.Each group in the polar plot is shown as one bin, and each polar plot shows α (i) in 36 angle bins.The length of each lobe in the histogram and its degree of darkness are proportional to the number of elements in α (i) (or number of pulses) that fall within a bin.On the edge of the round panel, we have a ring on which the degree of blackness reflects the total azimuthal distribution during the last 6 days of observations before the EQ.This "search light" indicates the most probable azimuths of ULF/ELF radiation and the corresponding plausible azimuths of a forthcoming EQ.They correspond to the darkest sectors of the ring.(1) The numerator contains the ratio of two horizontal spectral components P hh (NS component of magnetic field) and P dd (EW component).The denominator is the root mean square (rms) of the deviation of signal ellipticity.
The expression of β is given by equation ( 1).
( ) (2) Here P dh and P hd are cross-power spectral densities, and Im means imaginary part.Schekotov et al. (2007) have found an enhancement in the spectral ratio of P hh /P dd and a reduction in the standard deviation of ellipticity before an EQ, and the parameter introduced by equation ( 1) is proved to be most sensitive and reproducible to seismic shock.The ellipticity or the ratio of minor axis to major axis is defined by tan β.The sense of polarization is characterized by the sign of β; when β>0, the polarization is right-hand (RH), and β<0 means the left-hand (LH) polarization.The linear polarization is expressed by β=0 (Fowler et al., 1967).
The field power spectral densities, P hh , P dd and their cross-power spectral densities P hd , P dh are calculated by using Fourier transforms with frequency resolution of about 0.1 Hz.Spectral components in a frequency range from 0.1 to 30 Hz are taken into account.They are averaged over one-Hz intervals such as 0.1-1, 1.1-2, ……, 29.1-30 Hz, so that we have 30 spectral components in the present analysis.
A success in the application of this parameter ∆S is partly due to the fact that a majority of nearby EQs take place east of our station at Karymshiro, Kamchatka, Russia (Schekotov et al., 2007).In a more general case with the rotation of axes by some angle φ i , we can find a maximum of ∆S ( φ i ), and its radial component is directed to the source of radiation.See the details on the effect of axis rotation in ∆S ( φ i ), which have already been described in Ohta et al. (2013).
Thus the sequence of processing is as follows: 1) We divide the daily file into 24 1-hour intervals.
2) For every interval, we calculate spectral densities and find a maximal spectral component of ∆S with the rotation of axes.
3) We find the maximal value of ∆S among all previous values.
We simplified the processing by finding an optimal time interval.An example of results obtained by the algorithm described above is shown in Figure 5.The top panel indicates geomagnetic Kp index and the local seismicity index (Kls) which is given by Molchanov & Hayakawa (2008), where R is the epicentral distance (in km).
The evolution of the spectrum of ∆S for positions of axes when ∆S reaches a maximum value is shown on the 2nd panel, in which ∆S is calculated for each one-hour interval.Its daily maxima depicted on the 3rd panel by black bars and one-hour values are shown by gray bars.We can observe here a growth in ∆S on 8 March, just before the foreshock on 9 March and 3 days before the main shock.It was found that the main contribution to ∆S comes from an interval from 12 to 21 hours LT.Therefore ∆S was calculated only for this interval, and the result of this procedure is shown on the bottom panel.

Analysis Results
The dependence of ∆S on seismicity for the 5-month time interval from 1 January to 31 May 2011 is summarized in Figure 6.Seismicity is represented by the index Kls on the top panel and geomagnetic activity is expressed by Kp index on the 2nd panel.The figure exhibits a growth in ∆S just before the 2011 Tohoku EQ and a decline after the main shock.Then the evolution of spectrum of ∆S is shown on the 3rd panel and values of maximal spectral components are displayed on the bottom panel.
We here comment on some interferences in Figure 6.You notice horizontal line (constant frequency) signals at the frequencies of 8, 14, 20, 26 Hz etc., which are known as Schumann resonances due to the global lightning activity (Nickolaenko & Hayakawa, 2014).However, these Schumann resonances are considered to be interferences to the study of pulses in this paper, degrading the value of ΔS, so that we pay attention to the frequency range below the 1st resonance frequency (8Hz), between 1st and 2nd harmonics, between 2nd and 3rd harmonics, between 3rd and 4th harmonics etc. Figure 6 suggests that there is no correlation between Kp and ΔS during the large interval of our interest.Sunspots caused an increased Kp on 11 March, but it does not lead to any increase in ΔS.A solar proton event on 8 March is not a strong event with its proton flux=50, because we know that some effects are observable only for proton flux much larger than 1,000.So the growth in ΔS just before the EQ is likely to have nothing to do with geomagnetic and solar activities.
The same dependence and evolution of azimuthal directions of radiation during the last week before the 2011 Tohoku EQ are shown in Figure 7.It illustrates a sharp growth in ∆S on 8 March and a growth in intensity of radiation on 10 March.Moreover, the azimuth of this radiation is found to be directed to Japan, the region of epicenters of forthcoming EQs.Rectangular panels display the same as in Figure 6.Azimuthal distributions for seven days of observation including the day of EQ are shown on the bottom round panels.The corresponding azimuthal distribution for the day just before the EQ is shown on a separate large round panel located to the right of rectangular panels.
Figure 7.The illustration of a sharp growth in ∆S on 8 March and a growth in intensity of the radiation on 10 March.The azimuth of this radiation is directed to Japan, the region of epicenters of forthcoming EQs.The evolution of seismicity is represented by Kls index and that of geomagnetic activity is presented by Kp index on the top panel.Two bottom rectangular panels display the same as in Figure 6.Azimuthal distributions during seven days of observation are shown on the bottom round panels.The azimuthal distribution for the day just before the EQ is shown on a separate large round panel located to the right of rectangular panels

Summary and Discussion
This paper is considered to be a further extension of our recent paper by Ohta et al. (2013), in which we have found atmospheric ULF/ELF precursors by using the ULF/ELF observations at three stations in Japan for the 2011 Tohoku EQ.In order to compare with Figure 7 obtained from a remote station in China, we re-plot the corresponding results in Japan (Ohta et al., 2013) in Figure 8. Figure 8 illustrates the behavior of ∆S and azimuth distributions of the ULF/ELF signal observed by induction magnetometers at three field sites, Nakatsugawa (NAK) (geographic coordinates: 35.42°N, 137.55°E),Shinojima (SHI) (34.67°N, 137.01°E), and Izu (IZU) (34.64°N, 138.85°E) of the Chubu University network (Ohta et al., 2013).Top left rectangular panels show the temporal evolution of spectral maxima of ∆S (the same as in Figure 7) at those stations, and three rows of the round panels on the bottom of the figure indicate the evolution of the azimuth distributions of atmospheric ULF/ELF radiation.The top right panel shows the position of three observation sites (NAK, SHI, and IZU) and azimuthal ditsributions of ULF/ELF emissions on a particular day of 10 March, 2011 (one day before the main shock).
The color and dimension of a circle indicate the focal depth and magnitude of the EQ, respectively.The scales of depth and EQ magnitude are given by means of two subplots on the right-hand side of the map.
A comparison of this figure with Figure 7 enables us to come to the following major conclusions: (1) Precursors of ∆S are observed not only in Japan, but also at a remote station in China (with a large epicentral distance of about 4000km).
(3) Directions of the main lobe in azimuthal distributions of atmospheric ULF/ELF emissions are found to coincide approximately with the region of epicenters of the coming EQ.
(4) On the other hand, we can indicate some minor differences between the Japanese and Chinese results.
(5) The precursors in ∆S occurred on 6 March at all stations in Japan, but it was observed on 8 March in China.
(6) An additional day (6 March) exists in Japan when the main lobe of the azimuthal distribution is directed to the region of the forthcoming EQs at all stations.Unfortunately it is a little difficult to compare the results obtained in China and Japan, because the ULF/ELF systems in two countries are not exactly the same, but close to each other.Also, the observation was carried out under different measurement conditions.Due to strong interferences in Japan, we were forced to use the data only during the nighttime during a LT period from 1h to 5h.While in China, we had a possibility to choose an optimal time period (LT=12-21h) and frequency band of analysis because of lower interferences.These factors might be a possible reason of the above discrepany between Chinese and Japanese results.
The epicentral distance in this paper is about 4,000 km, but our previous study has detected a reliable growth of ΔS before an EQ (M=6.7) at a far distance of 915 km in Japan (Hayakawa et al., 2012b).One more point we have to mention here, is the previous work by Cohen & Marshall (2012), who have suggested no such atmospheric radiation in the higher VLF range (0.1-30 kHz) before the Tohoku EQ.But, their conclusion is not inconsistent with our present result, because their frequency range is obviously above ours.
An additional conclusion is that we can demonstrate that the use of ∆S initially proposed by Schekotov et al. (2007) and used subsequently by Hayakawa et al. (2012b) and Schekotov et al. (2008Schekotov et al. ( , 2013) ) is very effective in identifying atmospheric ULF/ELF radiation prior to an EQ, which will be of potential use in future in predicting an EQ.
Finally, we comment on the possible mechanism of seismo-atmospheric ULF/ELF emissions in the frequency range from units to tens Hz.Though it is quite uncertain at the moment why such low frequency radio emissions are generated prior to an EQ, Schekotov et al. (2007Schekotov et al. ( , 2013) ) and Molchanov & Hayakawa (2008) have suggested a few possibilities as generation mechanism of those ULF/ELF emissions as a seismo-atmospheric signature: the perturbations of the electric conductivity of the atmosphere caused by the ionized gas release (as suggested by Sorokin & Hayakawa (2013)) which may provoke not only usual lightning discharges, but also sprite-like-discharges between clouds and the ionosphere (e.g., Fullekrug et al., 2006;Nickolaenko et al., 2010).These discharges due to their substantially large dimensions may provide the most powerful radiation in the ULF/ELF range, and help us to explain the observation of ULF/ELF radiation at large distances as found in this paper.Moreover, the propagation loss in these frequency ranges is known to be quite small (e.g., Nickolaenko & Hayakawa, 2014).We can also suppose that this radiation may be caused by atmospheric gravity waves (AGWs) or infrasound turbulence excited by sporadic water/gas eruptions during a time interval of about 2 weeks around the EQ time (Molchanov & Hayakawa, 2008).
The last hypothesis may provide a possibility to explain 2-days delay of the peak of ∆S in China (if we assume that this delay is significant).This effect may be caused by discharges in the local atmosphere due to AGWs which came into this region in 2 days.It may occur in the vicinity of the point of observation in the region where there exist certain atmospheric conditions for thunderstorm activity (e.g.thunderstorm cloudiness).In our case its distance is small of the order of a few hundred kilometers, because the spectrum of ∆S has noticeable low frequency (<5 Hz) components particularly on 8 March (see Figure 5 and Figure 7).They should vanish for distant sources.Possible directions of the source are shown on one of round bottom panels (see Figure 7) which depicted the azimuth distribution of ULF radiation observed on 8 March.Most likely this source is located to the south.Let us estimate the horizontal velocity of AGW, V=4e6/(2*86400) ~ 23 m/sec, which roughly coincides with the propagation velocity of AGWs (Kuo et al., 2009).Pay attention to two bottom panels of Figure 5. ∆S visibly increases from 6 March (from 7 March in the bottom panel) to 8 March and it is decreasing during next days.AGWs (and the corresponding atmospheric disturbances) are approaching the region of our observation and then they are moving off.This possibility can be further investigated with the use of one more ELF station.
We have to mention that the mechanism of generation of ELF radiation before an EQ is poorly understood, so that it is needless to say that further extensive studies on the mechanism are highly required in future.

Figure 1 .
Figure 1.Frequency dependence of the magnitude of the ULF/ELF sensor transfer function

Figure 3 .
Figure 3. Map of the region of observation.The position of the magnetometer is shown by a black diamond and marked by a letter YS.EQs with magnitude (M) >6 which occurred during January-May 2011 are shown by color circles.Their color depends on depth and their dimension depends on the EQ magnitude.These correspondences for depth and magnitude (M) are shown by means of two subplots on the right bottom part of the map

Figure 4 .
Figure 4.An example showing the azimuth distribution of ULF/ELF emissions.The length of each lobe in the histogram and its degree of darkness are proportional to the intensity of radiation.On the edge of the round panel, we have a ring on which the degree of blackness reflects the total azimuthal distribution during the last 6 days of observations.It indicates the most probable azimuths of ULF/ELF radiation and the corresponding plausible azimuths of a forthcoming EQ.They correspond to the darkest sectors of the ring

Figure 5 .
Figure 5. Illustration of the processing algorithm.The top panel indicates geomagnetic Kp index and the local seismicity index Kls.The evolution of the spectrum of ∆S for positions of axes when ∆S reaches a maximum value is shown on the 2nd panel.Its daily maxima are depicted on the 3rd panel by black bars and one-hour

Figure 8 .
Figure 8.The evolutions of ULF/ELF field characteristics during the last week before the 2011 Tohoku EQ registered at three Japanese sites.Top left rectangular panels show the evolutions of ∆S at IZU, NAK, and SHI.Three rows of the round panels on the bottom part are the evolutions of the azimuth distribution at these sites.Top right panel is a map showing the position of sites, and azimuth distributions at each of the stations on 10 March, 2011