Seismogenic ULF/ELF Wave Phenomena: Recent Advances and Future Perspectives ()
1. Introduction
The new science field of seismo-electromagnetics has achieved remarkable progress in the last three decades [1] - [6] . Seismo-electromagnetics is defined as the study of electromagnetic phenomena associated with earthquakes (EQs) (but mainly before an EQ) for the final goal of short-term (lead time of about 1 week) EQ prediction, and it is agreed that these non-seismic electromagnetic effects do appear prior to an EQ. Various seismogenic phenomena have been found in the whole regions of lithosphere, atmosphere, and ionosphere, as observed from the ground and from space (satellites), but it is until recently, we come to a consensus that the ionosphere (not only the bottom, but also the upper F region) is extremely sensitive to pre-EQ lithospheric activity; see the book by Molchanov and Hayakawa (2008) [3] and reviews by Hayakawa (2011) [7] and Hayakawa et al. (2018) [8] for the bottom ionospheric perturbations as seen by subionospheric VLF (very low frequency)/LF (low frequency) propagation signals, and as for the upper ionospheric perturbations see Liu (2009) [9] and Liu et al. (2018) [10] as based on the observations by ionosondes and TEC (Total electron contents). Already the statistical correlations of those ionospheric perturbations both in the lower ionosphere and upper F region have been confirmed [9] [11] . Recent extensive studies on seismogenic ionospheric perturbations have been reviewed (e.g., [12] [13] [14] ). Furthermore, in addition to those ground-based observations, the subsequent observations by a satellite, DEMETER, have added more evidence on the seismogenic ionospheric disturbances (e.g., 6, [14] [15] [16] [17] ). As mentioned above, it seems that there have been established clear statistical correlations of ionospheric perturbations with EQs, though it is needless to say that the causality between the two has not been so identified, neither the mechanism of ionospheric perturbations (i.e., Lithosphere-atmosphere-ionosphere coupling (LAIC)) is poorly understood. So this seismogenic ionospheric perturbation is not dealt with in this paper. Another plausible candidate for EQ prediction seems to be DC geoelectric signal, the so-called SES (Seismic electric signal), which had the longest history and has been mainly investigated in Greece and Japan [18] [19] . However, this SES signal is likely to be a candidate for medium-term (lead time of about one month), so this topic is neither dealt with in this paper.
Except for the above two very promising physical items for short-term EQ predictors, there is another core group of seismo-electromagnetics, which are natural electromagnetic noises in different frequency ranges as possible candidates for short-term EQ predictors, which is the purpose of this review. Hayakawa et al. (2019) [20] have made an extensive review of those seismogenic ULF (ultra low frequency)/ELF (extremely low frequency)/VLF (very low frequency) electromagnetic noises including the following: 1) ULF/ELF transients (or Q bursts), 2) Anomalies in Schumann resonance (SR), 3) ULF magnetic field depression, 4) EQ effect on Pc 1 pulsations, 5) VLF electromagnetic emissions, 6) EQ effect on ELF/VLF whistlers, and 7) Lightning.
In our opinion, ground-based observations are likely to be much more favorable for realizing the real EQ prediction than satellite measurements, because we can trace the detailed temporal and spatial variations of any physical parameter only based on ground-based observations. Satellites can monitor mainly in-situ parameters within the ionosphere (i.e., ionospheric plasma, electric and magnetic fields, etc.) [15] [21] [22] [23] , but some of the satellites are aimed at monitoring the Earth’s surface parameters such as atmospheric temperature, brightness temperature, surface latent heat flux (SLHF), outgoing longwave radiation (OLR), thermal infrared spectral range (TIR), aerosols, etc. [6] [14] [24] [25] [26] . We do not deny that some of these parameters from space can be a plausible candidate for short-term EQ prediction. However, after the extensive consideration of several parameters in our previous review [20] , we have decided to pay the greatest attention to the following physical parameters with special reference to the estimation of their statistical correlation with EQs in terms of probability gain (PG) and their possibility to be used for the real short-term EQ prediction or scientific EQ prediction practice. In this updated version, we have selected the following four-wave topics.
1) ULF lithospheric radiation;
2) ULF horizontal magnetic field depression (as an indicator of lower ionospheric perturbation);
3) ULF/ELF atmospheric radiation;
4) SR anomalies (perturbation in the stratospheric and lower ionospheric perturbation).
The first topic was not treated in our previous review [20] , but we have added this in this review. These four parameters are highly likely to be important further candidates for real short-term EQ prediction in addition to ionospheric perturbations. One by one we will first review the essential points again in this version because all of these phenomena except the first one are not well recognized in the scientific society, and present the recent advances and future perspectives with special reference to their statistical significance in terms of PG in the practical EQ prediction.
2. Consideration of Real EQ Prediction: Estimation of Statistical Correlation with EQs in Terms of Probability Gain
When thinking of the real short-term EQ prediction of a future EQ (when (time), where (position) and how big (magnitude)?), the ground-based observation is much more useful and promising than satellite observations. The most important advantage of ground-based observation of any physical parameter is that it provides us with the detailed temporal evolution of a physical parameter. While, satellite observations of any physical parameter have been characterized by the disadvantage of discontinuous observation and the measurement depends on local time in a fixed area, so that it may be highly possible to miss some EQs or catch an incomplete ionospheric temporal-spatial evolution before some events. As for the ground-based observation, there are two different categories; that is, “local” and “integrated”. Each category has its own advantage and disadvantage. First, the local measurement is itself able to monitor an area just around the place of a sensor, so that this point is alternatively a disadvantage such that it is extremely difficult to accumulate the number of events. On the other hand, such an “integrated” measurement, the best example of which is the use of transmitter signals including VLF/LF (low frequency) transmitter signals or VHF (very high frequency) transmitter signals. For these transmitter signals, we can monitor any EQs taking place close to the great-circle path between a transmitter and a receiver, enabling us to accumulate the number of events much more easily than the local measurement.
Next problem will be closely related with the above issue on the number of events for any physical parameter: Statistical significance of anomalies of a physical parameter with EQs, such as PG or Molchan’s error diagram [3] [27] . In the real EQ prediction, it is useless to predict an EQ with magnitude of class 5, and our target will be EQs with magnitude greater than 6 or preferably greater than 6.5, which would lead to any significant damage to human lives and economic loss. Because [11] [28] [29] [30] offered some evidence on the scaling law of an effect size more pronounced for larger magnitude EQ events. Of course, when we decrease the EQ magnitude to increase the number of events as in scientific studies, it will inevitably deteriorate the statistical significance.
So, the main attention of this review is to pay particular attention to the future real EQ prediction, so that we will focus on the PG or so for the above-mentioned four physical parameters. Let us estimate the reliability of any phenomenon by using a conventional approach, which was developed to estimate the seismic precursor efficiency. We reproduce the definitions by Console (2001) [31] Molchanov and Hayakawa (2008) [3] :
• Target volume
is a volume in 3-D (three dimensional) space (time and 2 coordinates of the Earth’s surface) determined by the time of observation and geographical area of observation. Each EQ with preconditioned magnitude threshold or target event is depicted as a point in the volume
.
• Alarm volume
is a volume in which an EQ related to that precursor is expected.
• If an EQ occurs in the alarm volume, it is called a success (S).
• If an EQ occurs outside the alarm volume, it is a failure of prediction.
• An alarm that is not associated to any target EQ is called a false alarm.
If NS, NA and NE are the number of successes, the number of alarms and the total number of EQs in the target volume, then various performance metrics commonly considered in EQ prediction evaluation are as follows:
• Success rate = NS/NA is the rate at which precursors are followed by target events in the alarm volume.
• False alarm rate = 1 − NS/NA is the rate at which precursors are not followed by target events.
• Alarm rate = NS/NE is the rate at which target events are preceded by precursors.
• Failure rate = 1 − NS/NE.
• Probability gain
is the ratio between the rate at which target events occur in the alarm volume and the average rate at which target events occur over the whole target volume.
Generally, a precursor can be considered as reliable if it achieves a PG value greater than unity [31] . In our case of single station observation both the target volume and the alarm volume turn into the time intervals, hence, PG relation can be rewritten as follows:
(5)
where Te is the observation time and TA is the alarm interal. This PG is commonly used in the following estimation of the applicability of any physical parameter for EQ prediction, though recently some other parameters are used to evaluate the effectiveness of EQ prediction by any other metric, such as the receiver operating characteristics (ROC) technique by focusing on the area under the ROC curve (AUC) [32] .
3. ULF/ELF/VLF Studies
When identifying seismogenic ULF/ELF/VLF perturbations, the main obstacle are space weather effects (such as geomagnetic storms (or disturbances)) and meteorological effects. Electromagnetic fields (or waves) observed mainly on the ground lie in the ULF (f < 3 Hz), ELF (3 Hz < f < 3 kHz) and VLF (3 kHz < f < 30 kHz) bands. Depending on the field of research, different terminology is used to define the various frequency ranges, but we follow the frequency nomenclature employed in ionospheric and magnetospheric sciences. It is recently agreed that the observation of these ULF/ELF/VLF electromagnetic fields can provide us with important information not only on the upper atmosphere (ionosphere/magnetosphere), but also on the atmosphere and the Earth’s crust (lithosphere) (e.g., [33] - [38] ).
We first need to know the main sources of these ULF/ELF/VLF noises. In the field of atmospheric electricity, the most famous one is so-called SRs (e.g., [36] [37] [39] [40] ) covering the frequency range from 7 - 8 Hz to 30 - 35 Hz, which are known to be resonances in the Earth-ionosphere cavity, due to global lightning activity. The SR resonance frequencies are found to be very stable, so this SR can be treated as VLF/LF transmitter signals. The next source is ionospheric Alfvén resonator (IAR) which is the resonance in the lower magnetosphere and whose eigen-frequency lies in the range from 0.5 - 0.25 to 3 - 5 Hz, but its generation mechanism is poorly understood, either related to lightning (atmospheric electricity topic) or to magnetospheric effect (space physics) [34] [35] . The field-line-resonance and cavity mode eigen frequencies of the magnetosphere are below this frequency range since they cover the interval of 10−2 - 10−3 Hz which are even smaller than above resonance frequencies [35] . These ULF waves are definitely of magnetospheric origin including Pc’s (continuous pulsations) and irregular pulsations of Pi s due to the instability of ion Cyclotron waves in the magnetosphere [35] [41] [42] [43] . In the ELF/VLF range, there are conventional sferics (atmospherics) from lightning discharges and it is also known that there have been observed a few kinds of magnetospheric ELF/VLF emissions such as hiss and chorus, which are generated in the equatorial plane of the magnetosphere by plasma instabilities, propagated along the magnetic field line, to be observed on the ground [44] [45] [46] . Another popular ELF/VLF phenomenon is whistlers, which are originated from lightning discharges in the opposite hemisphere, propagate along the magnetic field line in the magnetosphere, penetrate through the ionosphere and observed on the ground in our hemisphere [47] [48] [49] [50] .
Though the main contribution to ULF/ELF/VLF ranges is from space physics and atmospheric electricity, it is found during the last three decades that electromagnetic phenomena with subtle intensity in these frequency ranges do take place mainly before the EQs (see, e.g., [1] - [6] ). In the following, we will focus on the above-mentioned four physical parameters in the ULF/ELF band to be identified from the background wave activity as described in this section.
4. ULF Lithospheric ULF Radiation
4.1. Historical Background
Since our previous review on seismogenic ULF/ELF/VLF wave phenomena by Hayakawa et al. (2019) [20] , there have been published several papers on ULF lithospheric radiation (e.g. [51] - [56] ). However, it is highly invaluable to describe historical events of the ULF magnetic field variations before huge EQs. The unprecedented ULF signature was first reported at Kodiak 2 h before the great Alska EQ (M (magnitude) = 7.2) of 27 March 1964 [57] . The following are the academic discoveries that local high-amplitude ULF signals were observed near the epicenters of violent Ms (magnitude (surface)) > 7 EQs at Loma Prieta, California [58] and Spitak, Georgia [59] [60] . The next one was observed for the 1993 Guam EQ (3 August, 1993, M = 8.0) [61] . The important indications from these EQ events can be summarized below [61] - [67] .
1) The main frequency range of lithospheric radiation is ULF around f = 0.01 Hz (period = 100 s), probably because of the skin effect of wave propagation in the lithosphere.
2) Magnetic field intensity exhibits an increase about one week (5 - 12 days) before an EQ, followed by a quiet period (quiescence), and a sharp intensity increase just before (less than one day) the EQ. The intensity ranges from a few nT to tens of nT.
3) The presence or absence of ULF anomaly is empirically expressed by the threshold of 0.025 R (epicentral distance (in km)) = M – 4.5 (see Figure 1); e.g., R = 70 km for M = 6, and R = 100 km for M = 7). The most fundamental point for this ULF radiation is that ULF radiation observation itself is a typical “local” measurement; i.e., a ULF observatory can sense only the region close to the observatory within 100 km or so.
4.2. Recent Advances
In recent years we know that several papers have been published to answer whether ULF lithospheric radiation can be considered as a reliable diagnostic for short-term EQ prediction [66] [67] [68] . By paying particular attention to some recent papers, we will try to answer this question.
First of all, we refer to important papers by Han et al. (2014) [69] . They have studied the ULF seismomagnetic radiation based on nighttime data observed at Kakioka Observatory (geographic coordinates: 36.23˚N, 140.12˚E) of Japan Meteorological Agency with sampling frequency of 1 Hz. Figure 2 illustrates their statistical results with superposed epoch analysis, which indicated that ULF magnetic anomalies are more likely to appear before sizable EQ events rather than after them; especially 6 - 15 days before an EQ within 100 km from the observatory. Then, Han et al. (2017) [70] evaluated their statistics with the use of Molchan’s error diagram shown in Figure 3, and they found that the PG of those ULF emissions is just around 1.6 against a Poisson model. Though this paper was criticized by Masci and Thomas (2015) [71] , their criticism was not so convincing.
Then, Warden et al. (2020) [72] have reproduced the results by Han et al. (2014) [69] , introducing an algorithm aiming to assess the significance of pre-EQ ULF anomalies. Also, they showed the sensitivity of their method to the choice of outlier rejection scheme, hypothetical depth, frequency band, and EQ catalog, and raised some concerns about the robustness of the results.
Next paper to discuss is the one by Yusof et al. (2021) [52] . They have used ULF data at 10 stations in the years of 2007-2010 in Southeast Asia, East Asia,
![]()
Figure 1. Summary of the seismo-ULF emissions in the form of earthquake magnitude (M) and epicentral distance (R). A white circle means the event with ULF anomaly, while a black circle, the event without ULF anomaly. The empirical threshold is indicated by a dotted line (0.025R = M − 4.5). After [62] .
![]()
Figure 2. Results of ULF geomagnetic anomalies. The blue, the black, and the red lines demonstrate the counts of EQ events, random_mean, and random_mean + 2σ, respectively. The gray bars and the green line show 5 day counts and corresponding random_mean + 2σ, respectively. Their values are given by the right vertical axis. Note that the count on the EQ day is not included in the calculation of the 5 day counts. The vertical black broken line (the 0 day) indicates the day when ES parameter is greater than 108. After [69] .
and South American regions, and finally 34 EQs with M greater than 5.0 have been analyzed. Extensive studies have been performed for these EQs, with the full use of polarization method [61] [73] and direction finding [74] [75] [76] . As the result, possible precursors of 20 EQs were identified with detection rate of
![]()
Figure 3. Molchan’s error diagram for EQ predictions using the geomagnetic anomalies as precursors. The leading time (Δ) is set as 11 days, and the alarm window (L) is set as 5 days. The black broken line gives the results of random guess; the black solid line presents the results of prediction based on ULF geomagnetic anomalies. After [69] .
about 60%. Figure 4 is taken from their paper, in which distributions of precursor presence are expressed in terms of a) epicentral distance (d, instead of R in Figure 1) versus M, and b) hypocentral depth (h) versus seismic index (KLS) [3] [75] (KLS = 100.75M/(d + 100)(d in km)). EQs are grouped into 3 groups; the bottom group below the 33rd percentile, 2) the middle group between 33rd and 66th percentile, and 3) the top group above the 66th percentile. The percentage of precursor presence, PP% is the highest for both M and d in their respective top groups, where PP% = 83% and 91% respectively. The result regarding M seems to be consistent with the previous assumption that greater M EQs (with higher KLS) emits higher intensity emissions, supporting a scaling law of effect size pronounced for larger M. While the high PP% for more distant EQs looks inconsistent with the summary (3) in the previous subsection, but this can be explained by the fact that most EQs having larger d coincidentally have higher M. From this figure, we can say that for an EQ with M greater than 6.7 there is an approximately 80% chance of having precursors.
We comment on the significance of these papers. Our personal feeling is that these ULF emissions must be associated with the generation of any electric current sources even though the mechanism is poorly understood. Taking into full account the extremely large wavelength in ULF, it is reasonable to consider that clear ULF emissions can be observable only for larger magnitude EQs, probably M greater than 6. Most of the sizable EQs treated in [69] had the magnitude less 6.0, and only about 10 EQs had the magnitude greater than 6. Even though it is useful to study EQs with M smaller than 6 scientifically and also as a scientific exercise of EQ prediction, it is too small a magnitude for us to make the real EQ
![]()
Figure 4. Distributions of precursor presences in terms of (a) epicentral distance, d against magnitude, M, and (b) hypocentral depth, h against KLS index. EQs were divided into 3 groups based on their percentile ranks, i.e., P<33, P33-66 and P>66. The percentage of precursor presences, PP% is shown for each group. After [52] .
prediction. So, it is easy for us to expect such a low PG of around 1.5, which is not sufficient enough for the real EQ prediction. At the time of writing their book, Molchanov and Hayakawa (2008) [3] already estimated the possible PG of these ULF emissions being just around 1, as shown as No. 4 in Table 1 (short-term precursors) and as No.3 in Table 2 (imminent precursors) for your reference (to be compared with different precursors).
4.3. Future Direction of ULF Emissions as a Reliable Diagnostic Tool
Because of the serious disadvantage of “local” measurement for lithospheric ULF emissions (though, of course, this is the most important advantage for short-term EQ prediction), it is extremely difficult to accumulate the number of ULF events, since we have to be located very close to an ULF observatory when an EQ happens. When trying to collect more events by decreasing EQ magnitude as in [69] , we suffer from an additional difficulty in including ambiguous ULF events with subtle intensity, leading to the insufficient PG of ~1.6 as in [69] . In order to get, at least, one step forward to an accurate EQ prediction with the use of ULF anomalies, we can suggest the following two criteria. Two criteria can be proposed for the acceptance of the observed ULF anomaly as an EQ precursor: 1) The reported existence of believable scientific testimony for ULF anomalies prior to, at least, some EQs, and 2) the existence of satisfactory physical models to describe the presence of ULF precursors.
In order to answer the first point (1), we think that the most suitable way to find more significant precursors to look for any ULF data close to the epicenter of big EQs with M greater than 6 (or preferably greater than 6.5 or so as suggested in [52] . However, we have to be very cautious about different factors. The prevailing arguments related to pre-EQ ULF emissions are mainly space weather effects (mainly geomagnetic disturbances). Few scientists [77] [78] [79] criticized the findings in some of the ULF studies as just a magnetic storm effect, but they failed to provide solid facts in support of their criticism. Still we have to watch carefully the relevant space weather conditions such as association with geomagnetic disturbances (or storms).
The epicentral distance of huge EQs should be close to a ULF observing station (preferably less than 100 km as indicated by the signal processing threshold), and this kind of ideal situations cannot be expected so often. We want to understand the whole view of wave characteristics such as wave polarization, wave spectra, wave intensity etc. for those ULF emissions seeming to be highly likely to be seismogenic. When the epicentral distance becomes larger like a few hundred kilometers, the ULF intensity is expected to be more subtle, so we are obliged to apply different kinds of sophisticated signal processing such as polarization, direction finding, etc. [62] .
Another useful suggestion is the establishment of net observations. We installed a ULF network with inter-station spacing of the order of 100 km such as done in the Kanto (Tokyo) district [62] during the Japanese Frontier Projects
![]()
Table 1. Short-term precursors (After [3] ).
*It is maximum time interval ahead of the EQ shock moment. **This estimation corresponds to the so-called alarm rate, i.e. the rate at which precursors are followed by large and nearby EQs (Ks > 1 or magnitude M > 6, distance D < 300 - 500 km). Such estimation is difficult to obtain. There are many reports on the subject, but many of them are a case study indeed (sign?). ***Abbreviations used here: WGM, water-gas migration through the ground medium: AGW, Atmospheric gravity waves; IT, ionospheric plasma-wave turbulence; EKE, electro-kinetic effect. Here PG is probability gain.
![]()
Table 2. Near-seismic precursors (After [3] ).
(though unfortunately this network already ceased), and recently as done by Chinese colleagues [55] . Beforehand we know highly EQ-prone regions based on the seismological studies (long-term prediction), so we are ready again to install such a ULF network in the relevant high seismic region. In normal conditions, we measure three field components of magnetic field, but additional measurements of electric field components are of special importance in estimating the whole wave characteristics of amplitude ratios and phase differences [67] [80] . This kind of multi-stationed network will enable us to perform sophisticated signal processing methods including principal component analysis [81] , singular value decomposition [82] [83] , direction finding (either goniometric or gradient [84] [85] or polarization ellipse [86] ), which will definitely increase the credibility of the observed seismogenic ULF emissions after eliminating the space weather effects.
All of the above-mentioned methods are expected to improve the credibility of seismogenic ULF emissions, and also for better understanding of their association with EQ activity, but no signal method can solve all of the above-mentioned problems and sophisticated combination of results by a few methods play an essential role in solving the problem.
An additional important suggestion is the use of AI (artificial intelligence) techniques such as machine-learning etc. in collecting more ULF events. Recently, those AI techniques have been proposed and applied to different parameters such as prediction of EQ location [87] , detection of seismic precursors in subionospheric radio data [88] . Further, there have been a lot of machine-learning applications to find ionospheric precursors from TEC data (e.g., [89] [90] [91] [92] ). In the field of geomagnetic data, Petrescu and Moldovan (2022) [93] have proposed a prospective neural network model for seismic precursory signal detection in geomagnetic field records. The algorithm does not take into account the physics behind the production of signals, but only picks up repeated features associated with the imposed labels. We need only a large enough dataset on which we train the machine-learning algorithm.
As for the generation mechanism (2), we have to confess that it is poorly understood, even though there have been proposed several hypotheses. Several seismogenic ULF emissions are believed to be originated from the local zones of EQs, because this is reinforced by the suggestions that mechanical deformation or microfracturing in the looming focal zones may give rise to pre- or co-seismic ULF emissions due to the following mechanisms: 1) electro-kinetic effects [94] [95] [96] [97] , 2) displacement of boundaries between high and low conductive crustal blocks [98] , 3) microfracturing electrification [99] , 4) crack propagation and the motion of charged edge dislocations [100] , 5) inductive effect resulting from the movement of a conductive medium in the Earth’s magnetic field [101] [102] . Further studies are highly required, in particular, on any attempts to associate the detected ULF emissions to EQ preparation process, or the generation mechanism.
Finally we emphasize the importance of critical analyses to any physical parameters including ULF emissions, even once there is obtained a statistical correlation with EQs. All of the above-mentioned analyses are based on the statistical analysis on the ULF amplitude increase with the use of standard deviation or so. However, we do not know whether the observed ULF emissions are really associated with an EQ, or an EQ precursor. This is a causality problem. There exist critical analyses in order to answer this question, such as fractal analysis (e.g., [103] [104] [105] ) or natural time analysis [106] . The use of critical analyses achieved a lot of success in identifying the association of EQ precursory anomalies with the EQs (for different parameters) (e.g., [106] - [116] ). Especially in the field of lithospheric ULF radiation, Hayakawa et al. (2015) [117] have used the natural time analysis to the ULF radiation before the 2011 Tohoku EQ and have confirmed the clear relationship of the observed ULF anomaly as the consequence of criticality in the lithosphere. This finding has been further confirmed for another EQ event [118] . So we recommend you to utilize this critical analysis to confirm the presence of seismogenic ULF emissions in future works.
5. ULF Magnetic Field Depression Phenomenon
This topic is not so popular even in our academic society, so we will repeat it here. The first report on the depression of ULF horizontal magnetic fields observed on the ground was reported by Molchanov et al. (2003, 2004) [119] [120] , and later confirmed on the basis of long-term observation [38] [75] [121] and some other case studies [122] [123] [124] . This phenomenon is characterized by the depression (or depletion) of the horizontal magnetic field component of ULF waves (mainly irregular pulsations (Pi’s) coming from the magnetosphere at night) (in the frequency (f
0.01 Hz) nearly the same as lithospheric ULF emissions) detected on the ground, which is likely to be attributed to the lower ionospheric perturbation [20] .
5.1. Initial Case Study
We studied our initial data from the observatory, Karymashino, Russia (geographic coordinates: 52.83˚N, 158.13˚E), which illustrates the temporal evolutions of the spectral density of magnetic field fluctuations during two months in the seven frequency bands from 0.003 Hz to 5.0 Hz. Z is the vertical component (Pzz) and G is the total horizontal magnetic field (G = Phh + Pdd). Though not shown as a figure, we could notice a correlation with magnetic activity Kp index, and a clear signature of correlation with Ks (seismic activity) (this KLS in the previous chapter) has not been seen in such a simple analysis.
In order to identify seismogenic ULF emissions from the ULF data, Hayakawa et al. (1996) [61] have proposed the use of polarization, which is the analysis of ratio of Z/G. The results for 7 months (from June 26, 2000 to February 25, 2001) are shown in Figure 5 taken from [20] . In the top panel we plot the geomagnetic activity (Kp) and seismic activity (Ks). Some correlation with Ks can be supposed at least for the frequency channels 2 and 3 and close to the date of large Ks values (these cases are indicated by vertical dotted lines). In order to check it we demonstrate all the cases in Figure 6, each case in the vicinity of the five most powerful EQs. Relying upon the results in Figure 5 and for simplicity we present only channel 2 (f = 0.01 - 0.03 Hz). It is obvious that nighttime values of Z/G exhibit an increase about 2 - 7 days before each EQ date.
An important question arises immediately: what does an increase of polarization mean? Is it due to an increase of Z component or a decrease of G component or both? We know that these changes practically synchronously. But because the value of Z is about an order smaller than G while they decrease together, Z quickly reaches a noise level of the sensor, or its reduction is limited by the interference, and then the polarization is likely to depend only on G. To confirm this assumption we present Z/G values and 1/G values after 1 day averaging in Figure 7 for the same case as in Figure 6.
Figure 8 illustrates the cross-correlation function (superposed epoch analysis) of the 2-day averaged local seismicity Ks and variation of ULF depression
. Here
![]()
Figure 5. Ks, Kp (top panel) and polarization Z/G in the different frequency ranges (subsequent 7 panels) during 7 months of observation at Karymshino station. After [20] .
![]()
Figure 6. Ks, Kp and polarization in the vicinity of the five most powerful EQs in the interval from 24/06/2000 to 25/02/2001. After [20] .
![]()
Figure 7. (a) The ratio of Z/G around the same EQs as in the previous figure after 1-day averaging. (b) The same as (a) but for 1/G values and in the corresponding scale. After [20] .
![]()
Figure 8. Cross-correlation function of the 2-day averaged local seismicity Ks and variation of the depression δ(1/G) at different frequency bands. After [20] .
and
is one month running mean value. This dependence shows that the depression maximum reliably exists before an EQ by about 2 - 4 days. This was the first result of a study of this phenomenon reported by Molchanov et al. (2003) [119] .
5.2. Statistical Result
This statistical study is based on Schekotov et al. (2006) [121] . Because we have accumulated the data for the 4-year period from June 2000 to June 2004, we were interested in whether we will see this phenomenon in the subsequent three years of observations. If so, it is crucial for the short-term EQ prediction whether this effect appears 2 - 4 days prior to an EQ. We were interested as well in the dependence of ULF depression on parameters of local seismic activity and the possibility of observing this phenomenon in other regions as well. At this stage we have defined a relative variation of depression δD, and the cross-correlation functions δD*Ks are calculated.
Because the effect is uncertain for EQ series, only isolated EQs are chosen for our analysis. We did not consider EQs with Ks < 1 and aftershocks in the interval less than one week after the main shock. As a consequence the number of isolated EQs or the first in a series with Ks ≥ 1 reduces to 39. One example is illustrated in Figure 9. Here is an example of variation of the field depression δD, Kp and Ks on the top-left panel and a map (top-right) with parameters of corresponding EQ which happened on 16/06/2003. Bottom panels are depression spectra for local noon (red) and midnight (blue). Maximum night-time depression occurs on the night of June 13, 3 days before the EQ. Especially the depressed intensity belongs to the frequency range of 0.02 - 0.05 Hz, while black line shows the inverse sensitivity of the magnetometer.
The cross-correlation function δD*Ks between Ks index and depression of the magnetic field δD is shown on the left panel of Figure 10, and the averaged variation of depression by the superposed epoch analysis is shown on the right panel. Both curves are found to have a clear maximum 3 days before an EQ.
In order to validate the statistical stationarity of the effect, the interval of observations is divided into 4 one-year intervals and the cross-correlation functions are calculated for each interval in the bottom panel of Figure 11. The clear maximum 3 days before the EQ is found for the 1, 3 and 4 year intervals as well as for the whole period of observations. As for the second interval, there take place additional maxima comparable with the “3-day” maximum. Thus, no meaningful effect is found for the second interval, but this is probably because of low level of the nighttime geomagnetic activity during the second interval (see the upper panel of Figure 11).
In [121] we had presented the dependence of depression on Ks, which indicated that Ks dependence becomes more or less clear for EQ magnitude value M > 5, and the correlation seems to be improved for EQ depth less than 50 km. Figure 12 illustrates our new figure on the statistical correlation of ULF depression with Ks index, based on the recent observations in the same Kamchatka observatory (number of events is 26) [125] . One circle corresponds to one ULF depression event, where the dimension and color of each circle indicate the EQ magnitude and depth, respectively. In the figure we use a linear interpolation of this dependence, and the red line shows the mean value of this relationship. While, three pairs of dashed lines indicate the boundaries of 70%, 80%, and 90% confidence intervals. This figure gives a further support to our previous results, and it clearly suggests a significant scaling law of effect size (ULF depression) more pronounced for larger Ks EQ events.
5.3. Further Examples
The ULF depression was observed also for the huge Japan EQ (M ~9.0) on March 11, 2011, which is a typical oceanic EQ of the plate type [124] . Figure 13
![]()
Figure 9. Top panels are an example of variation of the field depression δD, Kp and Ks and a map with parameters of the corresponding EQ. Bottom: depression spectra for local noon (red) and midnight (blue) for the seven days. After [20] .
![]()
Figure 10. Left: Cross-correlation between Ks of the selected EQs and field depression for the whole period of 4-year observation. Right: Variation of the field depression averaged over all the selected EQs by SPE method. After [20] .
![]()
Figure 11. Top: Ks index for EQs and year-mean field magnitude at the frequency range 0.03 - 0.5 Hz in the vicinity of ±0.5 hour around the local midnight. Bottom: cross-correlation functions for each year interval. After [20] .
![]()
Figure 12. Statistical relationship of ULF depression with KLS index. One circle indicates one ULF depression event, and the dimension and color mean the EQ magnitude and depth. The red line shows mean value of this relationship. Three pairs of dashed lines indicate the boundaries of 70%, 80% and 90% confidence intervals. After [125] .
![]()
Figure 13. ULF magnetic field depression for the 2011 Japan mega-EQ. The second panel refers to the station of MMB, the third panel, KAK, and the last panel, KNY. After [124] .
is their result of analysis of relative ULF magnetic field depression (δDep, exactly the same as δD in the previous subsection) at three Japanese ULF observatories (Memambetsu (MMB) (43.91˚N, 144.189˚E) in Hokkaido, Kakioka (KAK) (36.12˚N, 144.186˚E) in the main island, and Kanoya (KNY) (31.42˚N, 130.88˚E) in Kyushu). The definition of this relative depression is given in [123] [124] . The top panel indicates the temporal evolutions of Dst and the occurrence of EQs with M greater than 5. The figure shows that remarkable ULF depression is taking place on March 6, 5 days before the main shock at all stations, but the most significant depression is found to be observed at KAK in the main island, which is closest to the epicenter of the main shock, and so this anomaly is considered to be a precursor to this EQ.
Two additional huge EQs happened in the Kuril Islands (M~8) have been analyzed, and the quite similar ULF depression for these EQs has been confirmed [123] . Further studies have been carried out [38] [122] .
5.4. Recent Advances
Hayakawa et al. (2021) [126] have made the multi-parameters observations for the recent two successive huge EQ (with M around 7) happened in the off-shore of Tohoku area, probably as aftershock of the 2011 Tohoku mega-EQ. Even though the observatory of Kakioka is located about 300 km from the EQ epicenter, they have found extremely clear depletion of ULF horizontal magnetic field at f = 0.01 - 0.02 Hz, a few days before the first EQ in February and about one week before the second EQ in March. On the other hand, the presence of lithospheric ULF radiation in the previous section is very uncertain.
Then, Hayakawa et al. (2022) [127] have tried to find out whether the corresponding ULF depression effect for a much smaller magnitude (M = 5.9) EQ in Tokyo, to be compared with the above previous result for sufficiently large EQs (M~7). Figure 14 is taken from their paper, which presents the analysis results of ULF magnetic field at the same observatory at Kakioka about 80 km from the EQ epicenter. The frequency is 0.01 - 0.02 Hz. The top panel indicates the Kp and Ks indices. The second and third panels refer to the magnetic horizontal and vertical components respectively, while the fourth panel refers to the wave polarization as the ratio of vertical to horizontal component (Pz/h). And the bottom two panels indicate the ULF depression, and the last bottom one to the absolute value of ULF depression. Though it is very uncertain about the presence of significant lithospheric ULF emissions as in the 2nd to 4th panels, the ULF depression is an obvious precursor a few days before this sizable magnitude EQ. So, the covering area for ULF depression seems to be rather wider than that of ULF lithospheric radiation.
5.5. Summary of Observations
We come to the summary of basic properties of the effect of ULF magnetic field depression:
• A few days before an EQ in the vicinity of local midnight, there occurs a decrease in horizontal magnetic field fluctuations of nighttime irregular pulsations in possible association with EQs;
• It is especially noticeable at the ULF frequencies, especially 0.03 - 0.05 Hz;
• The absolute value of ULF depression seems to depend linearly on the local seismicity (or M), exhibiting a scaling law of more pronounced for larger M EQs;
• Despite of the same local measurement, the sensing range of ULF depression seems to be much larger than that of lithospheric ULF radiation.
Of course, it is needless to say that we need to accumulate the number of events by collecting more and more ULF data at a station close to the epicenter of huge (M ≥ 6.5) EQs, leading to better understanding including the physical process of this depression effect as will be discussed in the next subsection.
5.6. Consideration of the Mechanism of ULF Depression
We studied the fluctuations of geomagnetic fields in the frequency range of geomagnetic pulsations (0.01 - 0.1 Hz) at middle and low geomagnetic latitudes in highly seismic regions of Kamchatka (Russia) and Japan in this Section. The basic idea of this phenomenon is that the parameters of natural ULF geomagnetic noise can change significantly near the epicenter of a forthcoming EQ due to the LAIC process.
We have been treating nighttime ULF geomagnetic field variations, because the artificial (industrial) noise is minimum during nighttime, so we deal with nighttime irregular pulsations Pi’s (e.g. [35] [41] [42] [43] [126] [127] ). The fundamental basis for this phenomenon lies in the fact that the amplitude of
![]()
Figure 14. Bottom two panels refer to the characteristics of ULF depression (depression Deph and relative depression δDeph). Top panel refers to the magnetic (Kp) and local seismic activity (KLS). Avertical dotted red line is the occurrence time of our EQ. For the sake of comparison, we have reproduced the results of the lithospheric ULF radiation; the second and third panels referto Fh and Fz, and the fourth panel to polarization Pz/h. After [126] .
ULF geomagnetic fluctuations depends both on the amplitude of incident wave and the parameters during their propagation in the magnetosphere/ionosphere [128] . It is not reasonable to think about the EQ effect on the equatorial plane of the magnetosphere where those pulsations are mainly generated, so it is quite natural to attribute some changes in the ground-based ULF noise to the change in the ionospheric parameter of the propagation medium as in the case of subionospheric VLF/LF (low frequency) propagation anomalies, because we already know that the lower ionosphere is definitely perturbed before major EQs [129] - [137] .
Molchanov et al. (2004) [120] suggested that the depression of magnetic pulsations is caused by any ionospheric disturbances before violent EQs in the following ways. We shall consider two scenarios to explain the mid-latitude ULF geomagnetic noise depression:
1) The source is far above the ionosphere, i.e. in the magnetosphere. The ionospheric influence is through the variations of transmission coefficient of downgoing pulsations from the magnetosphere during the pre-EQ phase (Model I).
2) The variation of the wave vector (k) of the ionospheric fluctuations during an EQ preparation phase (Model II).
For the first scenario (Model I), incident “noise” (broadband) as hydromagnetic waves fall onto the ionosphere, which are registered on the Earth’s surface as natural ULF geomagnetic noise. The incident ULF wave from above is likely to be Alfvén waves (ion Cyclotron waves) [35] [128] . Following the mathematical formulations [128] [138] , the reduction in the transmission coefficient of the incident (downgoing) Alfven wave, or the ionospheric screening, grows with an increase in Pedersen conductivity or a decresase in Hall conductivity in the ionospheric dynamo region. So, on the reasonable assumption of the increase in Pedersen conductivity, Sorokin et al. (2004) [139] estimated the decrease in the transmission coefficient of incident Alfvén wave (or ground-based ULF wave intensity), which is found to be qualitatively consistent with the observational fact. This can be corresponded to a significant growth in electron density in the lower ionosphere by 5 - 10 times. This assumption is not so unrealistic, because the lower ionosphere is definitely found to be perturbed prior to an EQ with the use of subionospheric VLF/LF propagation anomalies [129] [132] [135] [136] . Much more quantitative estimations for varying ionospheric conditions are highly required, and it is now in progress. The origin of ULF depression is definite to be located in the lower ionosphere, whatever its mechanism is either a change in transmission coefficient of magnetospheric Alfvén wave (Model I) or the decrease of spatial scales of ionospheric sources (Model II).
However, the mechanism of why the ionosphere is perturbed before an EQ, is a challenge and is poorly understood, but we show the plausible hypotheses of this LAIC process so far proposed. Emergence of lower ionospheric perturbations prior to a large EQ has already been evidenced by means of subionospheric VLF/LF propagation anomalies [129] - [137] . Because these subionospheric VLF/LF signals are known to be reflected from the lowest ionosphere (D region at day and lower E layer at night), it is already believed that the lowest ionosphere is perturbed a few days to about a week before large EQs. Furthermore, the upper F region of the ionosphere is also known to be perturbed before the major EQs with the help of observational data on ionosondes, GPS TEC (total electron contents), satellite observations, etc. [2] [3] [6] [13] [14] . A lot of evidence has been presented with an enormous rate on the presence of perturbations in the ionosphere as presented in this paper with the use of GPS TEC observations because those GPS data are available openly worldwide [14] .
Hayakawa et al. (2004) [140] have already proposed a few possible hypotheses on the LAIC mechanism: 1) chemical channel, 2) atmospheric oscillation channel, and 3) electromagnetic channel, which are shown schematically in Figure 15. As for the first channel, emanation of radioactive radon induces the perturbation in the conductivity of the atmosphere, a change in the atmospheric electric field that leads to the ionospheric modification through the atmospheric
![]()
Figure 15. Schematic illustration of the lithosphere-atmosphere-ionosphere (LAI) coupling 1) chemical (+electric field) channel, 2) atmospheric oscillation channel, and 3) electromagnetic channel.
electric field [2] [140] [141] [142] [143] [144] . The second channel is based on the key role of atmospheric oscillations (acoustic wave (AW) or atmospheric gravity wave (AGW)) in the LAIC, and perturbations in the Earth’s surface (such as temperature, pressure, etc.) in a seismoactive region excite the atmospheric oscillations traveling up to the ionosphere and inducing ionospheric density perturbations [145] [146] [147] . The last mechanism, 4) electromagnetic channel, is that radio emissions (in any frequency range) generated in the lithosphere propagate up to the ionosphere and modify the ionosphere through heating and/or ionization. But this mechanism is found to be insufficient because of the weak intensity of lithospheric radio emissions [148] , and so the first and second mechanisms are likely plausible candidates for this coupling at the moment. Pulinets and others [2] [141] insisted on the first chemical channel as the most promising candidate for the ionospheric perturbations associated with EQs. That is, the emanation of radon was suggested as the main player of seismo-ionospheric perturbation, but there seems to be very little experimental (observational) evidence in support of their hypothesis, even though we know that there have been reports on the radon emanation itself (e.g., [149] ) as a precursor to an EQ. This chemical channel has been seriously criticized by Denisenko (2015) [150] , Prokhorov and Zolotov (2017) [151] , and Surkov et al. (2022) [152] because of its low efficiency. The most important signature associated with this first channel seems to be the generation of an electric field in the ionosphere in possible relationship with EQs. We give a brief explanation here. Chmyrev et al. (1986) [153] observed, using the Intercosmos-Bulgaria 1300 satellite, the anomalies in ULF (0.1 - 8 Hz) electromagnetic and quasi-static electric fields above the EQ epicentral zone and its conjugate region and also prior to an EQ. Such observations of seismogenic electric fields and associated phenomena within the ionosphere have been summarized by [2] [5] [142] . Unlike the generation of electric fields in the ionosphere, there have been no such significant changes in the electric field near the Earth’s surface. Both of these phenomena of conspicuous electric fields in the ionosphere and no noticeable changes in electric field near the Earth’s surface prior to an EQ would be an appropriate special subject to be discussed for the first channel.
Compared with the first chemical (+electric field) channel, a lot of experimental evidence has been recently accumulated in favour of the second channel (due to atmospheric oscillations) mainly by using subionospheric VLF data. Hayakawa et al. (2011) [147] have summarized different kinds of signatures of atmospheric oscillation effects in support of the second channel, such as the finding of many examples of enhancement of fluctuation in subionospheric VLF/LF data in the frequency range of AW and AGW before an EQ. In addition to the above indirect evidence, we have provided here the first direct evidence of AGW hypothesis using the VLF/LF Doppler observation [154] . We have no definite observational evidence on the correlation between information near the Earth’s surface and in the ionosphere, though we know very well that co-seismic ionospheric perturbations are apparently due to atmospheric oscillations (AW and AGW) (e.g., [155] [156] ) because the focal zone of an EQ induces a strain process near the Earth’s surface that would excite the AW and AGW propagating into the ionosphere. Korepanov et al. (2009) [157] made the first attempt in this direction by making full use of the data of ground-based pressure and magnetic fields and satellite-based plasma, but with using meteorological disturbances, to indicate the importance of AGW channel for the LAIC process. Nakamura et al. (2013) [158] have then extended the idea of [156] to pre-EQ effects, who have indicated some possible relationship between the Earth’s surface change and the ionosphere. Recently, Yang et al. (2019) [159] [160] have provided the first convincing evidence of the presence of AGWs in the atmosphere/stratosphere with the use of satellite data for a particular EQ (the 2016 Kumamoto EQ), followed by the similar finding for the 2011 Tohoku EQ [161] . Further evidence has been provided by Kundu et al. (2020) [162] and Politis et al. (2020) [163] . Especially, [163] has examined 12 huge (M ≥ 6.7) EQs, and found that all EQs except one supported the AGW channel, together with the critical analyses. Of course, it is necessary to perform further detailed studies in order to acquire observational facts in favour of the second channel.
An alternative hypothesis (electrostatic channel) has also been suggested by Freund (2009) [164] , who summarized his own results based on laboratory experiments. The discovery of positive holes charge carriers in crustal rocks, alongside electrons, opens a window of opportunity to study all these pre-EQ signals. He shows that various pre-EQ signals are the consequence of this single process: stress activation of electrons and positive holes in rocks. When the positive holes arrive at the Earth’s surface, they can cause a variety of effects including ionization of air at the ground-air interface, perturbations in the ionosphere, and distinct infrared (IR) emissions. When and how electric currents flow deep within the crust depends on the flow pattern. Under certain conditions, probably late in the EQ preparation process, the “battery circuit” can close. In this case, the electric currents flowing out of the stressed rock volumes can become very large, potentially on the order of millions of amperes, leading to powerful ELF/ULF emissions, but the penetration of electric field into the ionosphere meets the same difficulty as in the case of the first channel.
Further, Sorokin et al. (2020) [165] checked the previous theories and have proposed a new modeling of the influence of EQ preparation processes on the ionosphere through the electric field and electric current occurring in the global atmosphere-ionosphere electrical circuit. Their model is based on the generation of electric fields as a result of injection of charged aerosols into the atmosphere (i.e., electromotive force (EMF)), which could explain the generation of electric field within the ionosphere without any significant electric field near the Earth’s surface, though not yet confirmed.
When we have a sufficient intensity of electric field in the ionosphere, it would lead to different kinds of perturbations or instabilities [2] [142] , though not shown. Though the details are rather complicated, the generation of such an EMF might result in AGW instability in the ionosphere, the formation of field-aligned current and plasma irregularities, and many other related phenomena.
5.7. Future Perspectives
Even though no paper has been published on the PB estimation for this ULF depression effect, the results presented in this review may provide us with a lot of hope that the PB for this ULF magnetic field depression is rather large as shown in Table 1. So, such further statistical studies are required, and also more examples for huge (M ≥ 6.5, or preferably M ≥ 7) EQs should be collected.
6. ULF/ELF Transients
In the field of atmospheric electricity, a Q-burst (or an ELF transient) was first named in old days by Ogawa et al. (1967) [166] , which was known as an extremely intense impulsive ELF radiation, but its generation mechanism has been left unanswered for a long time as a mysterious event. In 1992 Williams et al. (1992) [167] gave a new insight to this Q-burst, as it is closely associated with unconventional positive cloud-ground lightning discharges leading to the generation of transient luminous events (TLEs) in the mesosphere [168] . In recent years ELF transients in association with lightning discharges and TLEs have been attracting a lot of attention in the fields of atmospheric electricity and space physics (e.g., [37] [169] [170] ).
In this section, we will present the quite similar kind of ULF/ELF bursts in possible association with EQs, in other words, we can call them seismogenic Q bursts (or impulsive ULF/ELF transients), which have been summarized in [38] [75] . The following descriptions are mainly based on the paper by Schekotov et al. (2007, 2013) [75] [171] and Hayakawa et al. (2019) [20] . The essential issue is how to identify seismogenic ULF/ELF transients from the conventional ULF/ELF sferics.
6.1. First Observations at Kamchatka (Case Study)
1) Observation and analysis
We use the ULF/ELF magnetic data observed at Kamchatka, so we explain the equipment there. Variations of the ULF/ELF magnetic field are being measured at the Karymshino observatory, Kamchatka, Russia since June 2000 with a 3-component induction magnetometer in the frequency band of 0.003 - 40 Hz, with noise level 0.16*f−1
(f, frequency) and conversion function 0.4*f V/nT in the frequency band f = 0.003 - 4 Hz and 1.6 V/nT in the band f = 4 - 40 Hz. The sensors for measuring the horizontal components of H and D are oriented along the magnetic meridian and transversally to it, and the Z means vertical. The parameters of all three sensors are identical with accuracy less than 3% in the conversion function and 2˚ in the phase characteristics. These discrepancies are corrected with the help of calibration circuits. The output signal is digitized with a sampling frequency of 150 Hz with the 24-bit data acquisition system (DAS) and written on a DAS hard disk. The data are copied and sent to the observatory for further analyses. See the details of the equipment in [75] [171] .
We first perform the preliminary routine data processing, which includes the correction of non-physical data and data gaps, filtration and decimation with the sampling frequency of 50 Hz. Then, to estimate the spectral and polarization characteristics of signals, we calculate their power spectral densities and cross-spectral densities: Phh, Pdd, Pzz, Phd and Pdh. In evaluating the average characteristics of signals we use their Fourier transforms with frequency resolution of ~0.2 Hz and time window of 30 min, but for impulsive transient signals we use a wavelet transform based on a 5-order complex Gaussian wavelet.
We compute the parameters of polarization ellipse using the conventional procedure [172] and define the angle
that the principal axis of polarization ellipse makes with the H-axis, or the orientation of polarization ellipse according to the following formula:
(1)
Re in the numerator means real part. Taking account of the signs of the numerator and denominator, the angle
is defined in the interval [−π/2, π/2] in the usual manner. We can also compute the parameters describing the ellipticity and sense of polarization ellipse in terms of the angle
.
(2)
Here Im means imaginary part. The ellipticity or the ratio of minor to major axis is defined by
and the sense of polarization by the sign of
, i.e.
> 0 or
< 0 as the polarization being right-hand (R) or left-handed (L) as measured when looking into the propagation wave. The ellipticity varies from −1 (L) to +1 (R) with zero for the linear polarization.
In the frequency range of ULF the amplitude of geomagnetic fluctuations decreases rapidly with frequency in such a way of ~3 orders lower at 1 Hz than at 0.01 Hz. At higher frequencies the level of natural geomagnetic emissions is much lower, but artificial interferences dominate at frequencies above 50 Hz. Thus it is reasonable to choose the frequency range from a few Hz to a few tens of Hz to look for weak signals associated with seismicity. Here are the data analyzed of more than 3-year monitoring of the magnetic field fluctuations in the frequency range 1 - 40 Hz in a seismically active region at Kamchatka. We compare different parameters of natural ULF/ELF emissions for the seismically active and quiet time intervals and select the field parameters, which are more sensitive to seismicity.
2) Analysis result
To select the best among the ULF/ELF parameters, the record was studied in the interval of about 1.5 month long around a seismic swarm in the middle of March, 2003 as a case study. We choose this interval to study a relation between the seismicity (Ks) and ULF/ELF field variations because of their remarkable behavior. Because the first half of the interval is seismically absolutely quiet, and the second one starts with the Ms = 5.9 shock on March, 15. This EQ is the first in the EQ series with slowly decreasing intensity. The second peak in seismic activity corresponds to Ms = 6 EQ registered on March, 19. Epicenters of almost all the EQs lie in the sea, east of the observation point. Seismic (Ks) and geomagnetic activity (Kp) and ULF/ELF field parameters are summarized in Figure 16.
It is seen from Figure 16 that we have an interval with the enhanced
ratio starting several days prior to an EQ and lasting several days afterwards. A similar but weaker effect is noticed in the power spectra of the field components. Such behavior of spectral parameters may correspond to a source located eastward (westward) from the observation point. Coincidence of the intervals with high seismicity and steep variations of the field parameters makes the assumption plausible on a physical relation between them. Namely, we can assume that an additional local source of ULF/ELF magnetic field fluctuations appears in the epicenter zone during the last stage of the EQ preparation and after the EQ. Characteristics of the ULF/ELF geomagnetic field in the frequency bands aside of SRs are shown in Figure 17. From top to bottom we plot the power spectral density of the horizontal components of the magnetic field (Phh, Pdd), their spectral ratio (
), ellipticity (tan(beta)) and ellipse orientation (theta). Figure 17 refers to the results for the frequency band 4 - 6 Hz below the fundamental harmonic of SR, while Figure 18 gives the same parameters for the frequency band 20 - 24 Hz above the third harmonic of SR.
In both frequency bands the spectral power of H-component and the ratio increase 3 - 4 days before the first shock of the EQ-swarm started on March, 15,
![]()
Figure 16. Upper panel indicates
(seismic index) and
(index of global geomagnetic activity). The second and third panels refer to the spectra of horizontal components of the magnetic field, H and D, and the bottom panel, the spectral ratio
. After [20] .
![]()
Figure 17. Characteristics of ULF/ELF noise in the frequency band 4 - 6 Hz for the same period of observation. From top to bottom: seismic index (Ks), index of the global geomagnetic activity (
) and power spectral densities, spectral ratio, the signal ellipticity, and the orientation of polarization ellipse. In the bottom two panels gray lines show the current value of parameters (tanβ, θ), medium black line is a running mean value (24 hours window), and upper and lower black lines show the ± 1σ (standard deviation) deviation from the mean value. After [20] .
![]()
Figure 18. The same as Figure 17, but for the frequency range of f = 20 - 24 Hz. After [20] .
2003. Simultaneously, we find decreases in the absolute value of ellipticity and the standard deviation (σ) of both ellipticity and ellipse orientation angle.
For the further analysis, we take the nighttime intervals of ±5 hours around the local midnight and calculate the mean values of the power spectral densities, their ratio and standard deviations (rms) of the ellipticity and ellipse orientation angle in the frequency band 4 - 6 Hz. Different combinations of spectral and polarization parameters were tested to select a parameter which is most sensitive to seismicity, which is a kind of tools to distinguish seismogenic effects from other noises. The behavior of different combinations of spectral and polarization parameters is summarized in Figure 19.
It is found in Figure 19 that sensitivities of all the parameters exceed that of
and are approximately equivalent. However, the lateral extension of a source and its explicit position influence the ellipse orientation rather than the ellipticity. Considering that the location and the size of a source can vary within a limited zone, we have chosen the following parameter
,
(4)
which presents the seismic influence better than any other parameters. The results of comparison of the efficiency of two parameters
and
are shown in Figure 20(a) for the interval of 2003.02.24-2003.04.06 and in Figure 20(b) for the interval of 2004.07.12-2004.08.08. Both parameters are found to increase with seismicity, but the time correspondence of
enhancements
![]()
Figure 19. Seismicity, geomagnetic activity and parameters of the magnetic field over nighttime intervals in the frequency band 4 - 6 Hz. The upper panel indicates the
seismic index,
(index of the global geomagnetic activity). Panels 2 to 5 refer to
,
,
,
. After [20] .
![]()
Figure 20. (a) Left: Seismicity and geomagnetic activity, middle: (Phh/Pdd − 1), and right: ∆S averaged over nighttime intervals in the frequency band 1 - 24 Hz for the temporal interval of February, 24 to April, 6, 2003 (time goes from bottom to top). (b) The same as (a), but for the interval of July, 12 to August, 8, 2004. After [20] .
with the groups of EQs is really amazing. The parameter
demonstrates several peaks in the seismically quiet intervals, which correspond to low amplitudes and/or unstable polarization of the signal. The relevant peaks are suppressed when using the parameter
. Thus, the enhanced
indicates the appearance of an additional signal with the polarization ellipse oriented along the magnetic meridian and the ellipticity stable at time scales of several hours. In a more general case of arbitrary direction to the source
and
should be changed by
and
, respectively, where indices n and t correspond to the directions perpendicular and parallel to the direction from the observation point to the source, respectively. The effect is recognized in the clearest way in the frequency band of 4 - 6 Hz.
6.2. Statistical Study at Kamchatka and a Trial of EQ Forecast with ULF/ELF Radiation
As a statistical study, the spectral density of the total horizontal power
for the whole period of observation is plotted in Figure 21. Two indices of local seismic
and global geomagnetic
activity are again shown in the upper panel, and
and
are given in the middle and bottom panels, respectively. They are calculated with 2-day averaging over night hours in the frequency band of 4 - 6 Hz. The total horizontal spectral power exhibits a typical seasonal variation with the maximum at local summer. On the other hand, the seasonal variation is not obvious in the variation of
, which suggests an evident correlation with seismic activity. Five intervals of high seismicity are clearly seen in the upper panel and each of them corresponds to the interval when the parameter
exhibited a significant enhancement. It is also important that the parameter
is not influenced by geomagnetic activity.
The influence of individual EQs on
is illustrated in Figure 21. It is seen from the figure that the EQs located to the east of the observatory at distances R < 300 km contribute mostly to the variation in
. Due to the specific distribution of EQs in the observation region, these eastward EQs are simultaneously closest to Karimshino as shown in Figure 22.
The analysis of data for individual EQs shows that the following EQ parameters correspond to noticeable changes in the magnetic field polarization and especially in
: depths < 50 km, magnitudes Ms > 5.5 (E > 1013 J) and epicenter distances R < 300 km. Finally, let us estimate the reliability of the effect by using a conventional confusion approach [31] as given in Section 2.
To estimate the time scale of the
pre-seismic variation and its threshold level
that maximizes the PG, we use a conventional superposed epoch analysis with the dates of main shocks with
> 1 taken as centers of time intervals. The analysis of the PG in dependence on the
shows that there occur noticeable changes in the interval
days around the EQ date and the optimal threshold value of
is
for chosen parameters of analysis. The results of superposed epoch analysis for 16 intervals with
> 1 major EQs and max(
) >
(where −15 <
< 0 and
is time in days from the EQ day) are shown in Figure 23.
Figure 23 shows that the effect of atmospheric ULF/ELF emissions (as a superposed epoch analysis) is almost symmetrical about the moment of the first
![]()
Figure 21. (Upper) Seismic and geomagnetic activity, (middle) ∆S and (bottom) the total horizontal spectral power for the whole period of observations. Field parameters are calculated in the frequency range of 4 - 6 Hz for local nighttime. After [20] .
![]()
Figure 22. Map of the main shocks in the vicinity of Karimshino (shown with a crossed circle and station code KRM). EQ epicenters are shown with circles. The color (in divisions of gray) corresponds to the maximal during the last 5 days before an EQ, and the size is proportional to the EQ magnitude. Only the EQs with > 1 and depths H < 50 km are shown. After [20] .
![]()
Figure 23.
variation ± 15 days in the vicinity of 16 main shocks with
>1 and max(
) >
. After [20] .
shock and the leading time is about 5 days. We calculate the daily averaged
and use the following rules:
• An interval with
exceeding the threshold level
is considered an alarm interval.
• The duration of the alarm interval is 5 days after the start day.
• The anomalies occurred from the main shock until 5 days after it are considered as being associated with the aftershock activity and are excluded from our consideration.
The application of these simple rules to the real observational data as a scientific EQ prediction practice (e.g., the year of 2001) is shown in Figure 24, in which black circles on the abscissa indicate EQs without any precursor, filled squares correspond to success and empty squares, false alarms. The results of the analysis of the effectiveness of the method for EQ prediction are summarized in Table 3 for each year.
The overall PG value for each interval of observation is presented in the last column of Table 3. All of them are found to be greater than unity with the average of PG = 2.68, which suggests the usefulness of this method as a short-term EQ precursor.
Not only the time but also the place of an EQ can be determined from the analysis of the characteristics of this radiation [173] . As an example, Figure 25 illustrates the result of data processing in anticipation of the start of the swarm on March 15. Here, in addition to determining the time of this start, we attempt to determine its azimuth. The three rectangular panels from top to bottom are: KS seismic index, spectrum of
and its averaged value
. With the round panels we show maps in azimuthal equidistant projection centered at the observation point. Here we also plot the histogram of azimuths of the radiation. In the bottom three circular panels are maps in azimuthal equidistant projection centered at the observatory for the 3, 2 and 1 day before the swarm and on the right top circular panel on the day it started. We find that on the day of March 13 when
reaches a maximum peak of the azimuthal distribution, it coincides with the azimuth of the EQ happened on March 15.
Let us summarize basic properties of the considered ULF/ELF radiation:
• A few days before an EQ there appear the impulsive emissions in a frequency range from a few to a few tens of Hertz;
• The combined characteristic of the field
proved to be most sensitive to this radiation;
• Lead time in
depends on the frequency and it varies from 1 - 5 days in units of Hertz up to 1 - 2 weeks for tens of Hertz;
• On average, it reach a maximum three days before the EQ;
• Azimuth of the radiation source calculated for impulsive component of the magnetic field is approximately equal to the direction to the epicenter of future EQs;
• A statistical relation with EQs demonstrates the usefulness of the method as short-term EQ precursor.
![]()
Figure 24. Trial of the EQ prediction technique for an interval from January 22, 2001 to December 30, 2001. Upper panel: EQs with
> 1 (hexagonal stars),
daily indices, lower panel:
. Alarm intervals and EQs are shown by markers above the horizontal axis. Black squares correspond to success, and empty squares correspond to false alarms, and black circles indicate EQs without precursors. After [20] .
![]()
Table 3. Summary of
efficiency as an EQ precursor for 3 years.
![]()
Figure 25. Three rectangular panels from top to bottom are KS seismic index, spectrum of
and its averaged value
. With the round panels we show maps in azimuthal equidistant projection centered at the observation point. The histogram of azimuths of detected radiation is also shown. The three circle bottom panels illustrate the histogram of azimuthal distributions for 3, 2, and 1 day before the first EQ and in the top right panel on the day of the EQ. After [20] .
6.3. An Example for the 2011 Tohoku EQ
We have presented a statistical study on the correlation between linearly polarized ULF/ELF bursts with EQs on the basis of Kamchatka observations, and here we will present our latest result of those ULF/ELF bursts for a particular recent and disastrous 2011 Tohoku EQ, which was based on the paper by Ohta et al. (2013) [174] .
1) Observation system and observation network
The geomagnetic data used here have been obtained by the Chubu University ULF/ELF network which consists of three observatories; Shinojima (abbreviated as SHI; geographic coordinates, 34.67˚N, 137.01˚E), Nakatsugawa (NAK, 35.42˚N, 137.55˚E) and Izu (IZU, 34.64˚N, 138.85˚E) [175] . Figure 26 illustrates the relative locations of three ULF/ELF observatories and the epicenter of the 2011 Japan EQ (11/03/11, the biggest orange circle). Also, we have plotted, in Figure 26, one representative observatory of Kakioka (KAK), and some other foreshocks (09/03/11) and aftershocks (11/03/11).
At each observatory we measure the magnetic field changes (H, D and Z components) by means of three orthogonal magnetometers in the frequency range of 0.1 - 24 Hz. The magnetometer is an induction coil sensor (similar to those used in Kamchatka), and the receiver attained a high sensitivity of about 0.05
at the frequency of 10 Hz. The details of the equipment are described in [174] [175] . The data observed at each observatory are regularly sent to the master station of Chubu University (at Kasugai).
Horizontal components of magnetic field are digitized at the sampling frequency of 100 Hz with the use of the 16-bit DAS and those data are stored on a hard disk. Those data are transmitted to the master station through a telephone line or an internet.
Here we describe the procedure of ULF/ELF magnetic field analysis in order to detect any seismo-atmospheric electromagnetic radiation and to determine the azimuth of its source, though some parts are repetitions of the previous section.
However, the preliminary routine data processing was applied before the main analysis. That is, this process includes substituting the interpolated data for short (several points) data gaps leading to some errors in DAS, band-pass filtration by means of 4-order Butterworth filter with cut-off frequencies 0.1 and 24 Hz. Two-directional filtration was applied to prevent a time shift of data.
We take exactly the same analysis procedure used in the previous section dealing with the Kamchatka data.
2) Observational results
Figure 27 illustrates the result of the spectrum ∆S(f) at a particular station of NAK during a period from March 4 to March 9, 2011, covering the date of a huge foreshock on March 9. First, we look at the results in the top rectangular panel. The top panel indicates the local seismicity index (KS).
In the top panel of Figure 27, there is one day of higher KS: a foreshock on March 9 (Ms = 7.3) (09/03/11 in Figure 26). In the analysis we have used only the local nighttime data in the JST (LT) from 0.5 to 5.0 h (total duration of 4.5 hours) where we expect the minimum local noise. The middle rectangular panel
![]()
Figure 26. Relative location of our ULF/ELF observing stations (Nakatsugawa (NAK), Shinojima (SHI) and Izu (IZU)) and the epicenter of the 2011 Japan EQ (11/03/11) (biggest circle). For the sake of comparison, one geomagnetic station at Kakioka (KAK) and some other foreshock (09/03/11) and aftershocks (two circles with 11/03/11). After [174] .
![]()
Figure 27. Detection of ULF/ELF radiation at NAK and determination of azimuth of the source before the 2011 Japan EQ. The top plot of the top left rectangular plot illustrates the relationship of local seismicity index (KS) with ∆S(f) and ∆S(fs). The second panel illustrates the temporal evolution of frequency spectrum of ∆S(f) (more black for more intensity) and the bottom panel indicates the temporal evolution of ∆S (fs = 9 - 10 Hz) during the period of March 4 to March 9. The top right panel illustrates the overall azimuth plots during the period of March 5-March 9, together with EQs with M > 5. The bottom three panels illustrate the polar plot of azimuth distribution of ULF/ELF radiation on March 6, March 7 and March 8, respectively (from left to right). After [174] .
illustrates the temporal evolution of frequency spectrum of ∆S(f), in which stronger intensity is indicated with darker black. We can find that on March 6 there is a remarkable enhancement of ∆S at the frequency of 9 - 10 Hz. This result is illustrated on the bottom rectangular panel where we show the temporal evolution of ∆S in this frequency range.
Here we comment on the effect of magnetic storms on the behavior of ∆S. Although this point was already discussed in our previous papers [124] , we repeat only the essential point here. Though the temporal evolution of Dst index as a measure of geomagnetic activity is not shown here, we know that the time of a minor magnetic storm is not coincident with our peak in ∆S and also the value of Dst is close to null on March 6 when ∆S is maximal [124] . So our peak of ∆S is very likely to be seismogenic.
Figure 28 is the comparison of ∆S at the three observatories of NAK, SHI, and IZU during a much longer period of February 1 to March 14, 2011. It is found from this figure that ∆S at all three observatories exhibits sharp maxima on the same day of March 6, which is 3 days before the March 9 first strong foreshock and 5 days before the March 11 huge EQ. The peak at NAK is conspicuously enhanced because of lower electromagnetic noise there, and the sharp peak on March 6 is still very remarkable at SHI. The electromagnetic noise environment at the third station of IZU is not so good enough that we expect a lot of fluctuations in the variation of ∆S. As the conclusion, ULF/ELF radiation detected by means of ∆S appeared 3 days before the first foreshock, which is indicative of the beginning of seismic activity. This is in agreement with the conclusions of our previous statistical studies and gives us a possibility to estimate the time of a forthcoming EQ.
Next we determine azimuthal distributions of the radiation. An example of their presentation is shown on the round panels of Figure 27. The distribution of α is represented by an angle histogram, which is a polar plot (as in the bottom circular panels of Figure 27) presenting the distribution of α values. Each group in each polar plot is shown as one bin, and each polar plot shows α 1) in 36 angle bins. The length of each lobe in the histogram and its degree of darkness is proportional to the number of elements in α 2) that fall within a bin. Examples of daily plots on three days (March 6, March 7 and March 8) are given in the bottom three panels of Figure 27. A summary of azimuthal distributions of ULF/ELF radiation for the last 5 days of observations is shown by a degree of blackness on the ring which is placed on the top right panel of Figure 27. Its most dark sectors are found to coincide roughly with the azimuths of probable forthcoming EQs. Their limits or probable errors are shown by dashed lines which cross the point of observation.
The azimuthal distribution of ULF/ELF radiation recorded on March 6 when the seismogenic ULF/ELF is strongest is shown on the map of our interest in Figure 29. The size of lobes and degree of blackness as well as in Figure 27 are proportional to the pulse flux density of ULF/ELF radiation. There are illustrated positions of observatories (NAK, SHI, IZU) and EQs with M > 7 occurred from
![]()
Figure 28. Temporal evolution of ∆S (9 - 10 Hz) at three observatories of NAK, SHI and IZU during a longer period of February 1 to March 14. After [174] .
![]()
Figure 29. Azimuthal distributions of ULF/ELF radiation recorded on March 6 at three observatories of NAK, SHI and IZU. Positions of observatories are shown by diamonds. Magnitudes and depths of EQs are represented by size and color of circles. After [174] .
March 6 to March 11, and the magnitude and depth of EQs are represented by size and color of circles. Despite of a long scattering of seismic disturbance, the maximum of azimuthal distributions are found to be roughly directed to the epicenter of the EQ.
By using the ULF/ELF data in a frequency range of 0.1 - 24 Hz, we have examined whether there is observed any ULF/ELF precursor to the 2011 March 11 EQ. The following facts have emerged from the present analysis, which is a kind of reconfirmation of the previous summary in Section 5.2.
1) The combined characteristic of the magnetic field, ∆S is again proved to be extremely useful in figuring out seismo-atmospheric ULF/ELF radiation.
2) The temporal evolution of ∆S is found to be peaked on March 6, which happened 3 days before the first foreshock which is indicative of the beginning of seismic activity.
3) The frequency of the maximum ∆S(f) is observed in the vicinity of the first SR.
4) The azimuthal distribution of ULF/ELF radio emission is found to be coincident approximately with the position of the main shock region.
Consequently, we come to the conclusion that the ULF/ELF radio emissions are highly likely to be generated as a precursor to the 2011 March 11 huge EQ.
6.4. Recent Advances
Prospective short-term EQ prediction has been continued at Kamchatka [125] [176] with the main use of this ULF/ELF electromagnetic radiation together with the ULF depression as discussed in the previous Section 4. The three parameters should be determined in the short-term EQ prediction: 1) when (time), 2) where (position) a next EQ is coming with 3) how big (magnitude). The source location (or azimuth because of a single station) of ULF/ELF atmospheric radiation gives us an approximate estimate of the EQ epicenter. Then, we estimate the local M in consequence of its statistical dependence of ULF depression on EQ magnitude (as summarized in the previous section). The date of a coming EQ is the easiest parameter to determine by the statistical dependence of delays of EQs relative to the dates of their precursors. The result of application of this method to real magnetic data is illustrated by the official prediction processes during a period of March-May 2016, and further limits and possible errors of the method as well as methods to enhance the reliability of this prediction are discussed. In order to solve the disadvantage of a single-stationed observation in Kamchatka, Schekotov et al. (2021) [177] made a further proposal to estimate the position of a coming EQ. They have studied the spatial statistics of local EQs in relation to Kuril-Kamchatka and Aleutian trenches. It follows from this statistics that more than 90% of events with M of more than 5 occur in the gap ~150 km west of the Kuril-Kamchatka trench and north-east of the Aleutian trench. Combining this with the direction of observed atmospheric ULF/ELF radiation, they have obtained the location of EQ epicenter of a next EQ.
Hayakawa et al. (2021) [126] have made the multi-parameters observations for the two successive huge (M~7) EQs in the offshore of Japan in February and March 2021. Even though the importance of multi-parameters observation has been strongly emphasized in [6] , very few papers have been published based on multi-parameter observation covering the lithosphere, atmosphere and ionosphere because it is extremely difficult to perform this interdisciplinary approach. Among eight physical parameters examined in our previous review [20] (lithospheric ULF radiation, ULF depression, atmospheric ULF/ELF radiation, stratospheric AGW (atmospheric gravity wave) in [20] , lower ionospheric perturbation (with subionospheric VLF/LF propagation anomalies), upper ionospheric (F region) (foF2 at Kokubunji), TEC observation), we have found that there are two obvious precursors; one is ULF depression at Kakioka (epicentral distance ~70 km), and second is atmospheric ULF/ELF radiation as observed at Nakatsugawa by Chubu University ELF network. Another important paper has been published by Hayakawa et al. (2022) [127] , again based on multi-parameter observations, with a target of a much smaller (M = 5.9) EQ in Tokyo on 7 October, 2021, in which a new parameter has been investigated, i.e. meteorological data (temperature, and humidity) at many stations from Japan Meteorological Agency. As the result of analyses, being the same result as the above paper, two obvious precursors have been identified; one is ULF depression and the second is atmospheric ULF/ELF radiation. These studies have further provided us with further evidence on the potential usefulness of this atmospheric ULF/ELF radiation in the future studies of EQ prediction.
6.5. Modeling of Seismogenic ULF/ELF Emissions
The generation mechanism of seismogenic atmospheric ULF/ELF transients is most poorly understood. We have discussed two possibilities closely related with LAIC or exactly the same as LAIC: 1) perturbations of the electric field induced by the pre-seismic ionized gas release (e.g. radon), and 2) AGWs or the infrasound turbulence excited by sporadic water/gas eruptions or by foreshocks and aftershocks during a time interval about 2 weeks around the EQ time. Based on the discussion on seismo-ionospheric perturbations (Subsection 4.6.), we think that the first possibility is less probable than the second, but the first process is likely to cause discharges in the atmosphere, which seems to be consistent with the impulsive nature of our observed ULF/ELF radiation. So, more extensive works are highly required in future. Recently Schekotov et al. (2022) [177] have examined the correlation of anomalies in atmospheric parameters (temperature and humidity) (closely related with the first process) and ULF/ELF radiation, and they have come to a conclusion that radon emanation in the chemical channel is not the origin of atmospheric ULF/ELF radiation.
Mareev et al. (2002) [178] and Molchanov et al. (2004) [119] considered second possibility and showed that for the typical periods 0.5 - 1 hours and horizontal scales of the disturbances of about λ~1 km AGWs propagate in the direction about 5˚ - 7˚ to the horizon. For distances 100 - 200 km corresponding altitudes are 10 - 20 km, and the resulting density fluctuations are about 1%. But the principal question about the mechanism of generating observed phenomena due to low amplitude disturbances is still unclear.
An alternative mechanism is also associated with the enhanced plasma concentration above the EQ preparation zone, but this was found to be not so attractive [20] .
6.6. Future Perspectives
Since the generation mechanism of this ULF/ELF radiation is most poorly elucidated, it needs further extensive works, even though the statistical correlation has already been established with rather high PG as shown in this paper and also in the book by Molchanov and Hayakawa (2008) [3] (No. 8 in Table 1), which indicate a potentiality of this atmospheric ULF/ELF radiation as short-term EQ predictor. Unlike the local measurement of ULF lithospheric radiation (section 3), this ULF/ELF transients can propagate over longer distances, so that the covering range for this phenomenon is a few hundred km up to a thousand km. The two parameters (when and where) of prediction of a coming EQ are rather easy to estimate, but the dependence of any parameter of this phenomenon (intensity, ΔS etc.) on EQ magnitude should be investigated extensively.
7. Anomalies in SRs
We review in this section the impact of seismic activity on the SR signals. In spite of the fact that the phenomenon of the global electromagnetic resonance (SR) is one of the most popular subjects in the atmosphere electricity, this area of research has the shortest history among the issues addressed in the present review. The unusual SR signals were noted for the first time in the records of the NAK observatory in connection with the powerful EQ in Taiwan so-called Chi-Chi EQ (Hayakawa et al., 2005 [179] ). The subsequent six-year statistical investigation of the NAK records related to the shocks of the magnitude exceeding 5.0 has confirmed the existence of the EQ-SR link [180] . The similar effect of intensifying higher SR modes was detected at the Moshiri observatory in Hokkaido, and it was associated with the Ping-tong EQ, which also took place in Taiwan [181] . In 2009, Ohta et al. (2009) [182] reported a new type of anomalous SR signals for two EQs that occurred in Japan.
Evident association with the EQs of anomalous signals observed in the SR band has recently attracted attention of many researches in the world [183] [184] [185] [186] . We expect that further observations will provide us with many novel experimental facts.
7.1. First Evidence of Anomalies in SR Phenomenon Associated with EQs (Case Study)
The NAK observatory monitors three components of ELF magnetic field (Bx, By, and Bz) using the induction magnetometer with three mutually orthogonal coils (see Section 5.3). The field component Bx is oriented along the north-south bearing, the By component corresponds to the east-west direction, and the Bz is the vertical field. The spectra of these fields are obtained using FFT procedure applied to 1024 samples. This corresponds to the time segments of 10.24 s or to the frequency resolution of 0.097 Hz (~0.1 Hz). We concentrate on the amplitude and the phase of two orthogonal horizontal magnetic field components Bx and By.
The following SR anomaly was discovered initially in Japan described by Hayakawa et al. (2005) [179] in association with a huge Taiwan Chi-Chi EQ (M = 7.3), which occurred on Sep. 21, 1999. The spectrogram in Figure 30 depicts the dynamics of ELF waves at the NAK observatory around the time of the EQ. Temporal evolution in detail of the spectrum power density is shown in Figure 31 where the specific narrow frequency band is shown ranging from 26.36 to 26.56 Hz.
The following notable results were formulated in this case study:
1) The fourth SR mode was strongly intensified in comparison with the regular conditions. The power spectral density usually decreases with the SR mode number n (the first mode peak is the highest one, and the intensity declines with n). The anomaly is a clear elevation of the spectrum in the vicinity of the fourth SR mode. This anomaly emerges about a week prior to the EQ main shock.
2) The substantial frequency shift of about 1 Hz is present against the usual fourth peak position in the By field component of magnetic field. This component is sensitive to radio waves arriving in the meridian plane or along the south-north direction coincident with the direction toward the future EQ.
7.2. Statistical Study of NAK Data (Six-Year Records of EQs in Taiwan)
The paper by Ohta et al. (2006) [180] treated the SR anomalies at NAK in the plus/minus two weeks in the vicinity of a EQ in the Taiwan area, provided that the magnitude was M > 5.0. The ELF records at the NAK observatory covered the six-year time interval from 1999 to 2004. Position of the main shock for every EQ is shown in Figure 32 by the circles of two kinds: the red circles with white character and the white circles with the black characters. The first set of circles corresponds to the EQs associated with the SR anomalies recorded at the NAK observatory in the Japan. The second group marks the EQs passed unnoticed in the SR records. The following results might be outlined. During the six-year observations, the 33 EQs took place in Taiwan and its neighborhood, and all of them had the magnitude exceeding 5.0. Four of these EQs passed with no SR disturbances at the NAK observatory. Thus, we have 29 EQs (out of 33) with the SR anomalies. Seven EQs associated with the anomalies have occurred inland (the fault type). The rest 22 EQs have struck in the sea (oceanic type). The following conclusions were made from these statistical data.
1) The goal of the data analysis was in demonstrating the obvious link between the unusual SR signals detected regularly in Japan at the NAK observatory and the EQs occurring in Taiwan.
![]()
Figure 30. Dynamic spectrum of ELF waves at NAK in the frequency range up to 50 Hz just around the Chi-chi EQ in Taiwan. The EQ (M 7.6 and depth (D) of 30 km) happened at 01 h47 m LT on 21 September. We can easily identify an anomaly in SR by comparing the spectra before and after 16 September. The SR fourth harmonic is extremely enhanced before and after the EQ. After [20] .
![]()
Figure 31. Temporal evolution of the intensity in the frequency range from 26.36 to 26.56 Hz (about 1Hz shifted from the conventional SR fourth harmonic) in relation to the main Chi-chi EQ and subsequent EQs (all aftershocks). After [20] .
![]()
Figure 32. Location of EQs (with M greater than 5.0) in and around Taiwan. Red circle with white character (EQ event number) indicates that this EQ is accompanied by the SR anomaly at NAK. White circle with black character, refers to the event without SR anomaly at NAK. After [180] .
2) The anomalous signals in the SR band were observed for all seven EQs that took place inland the Taiwan.
3) The SR anomaly was found only twice among the 22 oceanic EQs. These two events (No.9 and No.20 in Figure 32) were specially addressed in this case study. It appeared that No.9 was the shallowest among the 22 events. The EQ No. 20 had the highest magnitude among the 22 others.
Thus, we arrive to the conclusion that there is an actual link between the SR anomalies observed in Japan at the NAK observatory and the EQs in Taiwan. The disturbances in the SR spectra are always present in the NAK records when an inland powerful EQ takes place in Taiwan. Some underwater EQs in Taiwan are able to cause the SR anomalies, provided that either their magnitude is great or their depth is small.
7.3. Modeling of Anomalies in SR Using Solution of Scattering Problem
We suggest an explanation below for the observational results in Japan mentioned in the previous sub-sections. The explanation implies the wave scattering and wave interference.
Two characteristic heights (h1 and h2) are introduced [186] for description the ELF radio propagation including the ULF/ELF transients of Section 5 and the SR anomalies in Section 6. The conductivity currents become equal to the displacement currents of a given frequency at the height h1, therefore, it is called the electric height. Physically, the electric height is the altitude where the atmosphere transforms from a predominantly non-conducting medium into the conducting one. The amplitude of electric field in the incident radio wave starts rapidly decrease with altitude above the height h1.
The height h2 is regarded as magnetic height. Here, the other physical process takes place: the radio wave proliferation changes from the propagation (described by the wave equation) to the spreading by diffusion (the heat equation). Only the magnetic field is able to reach the height h2.
In the SR frequency band, the characteristic heights are approximately equal to h1 ~55 km, and h2 ~95 km [36] [37] [187] . One may observe that the altitude interval of the ELF radio wave reflection lies below the reflection height of VLF waves. Thus, the ionosphere region significant for the ELF propagation might be disturbed to a greater extent by the seismic activity then the interval important for VLF waves.
The horizontal size of the disturbance is found from the equation
where radius R is measured in km and M is EQ magnitude [188] . One obtains R » 2000 km for the Chi-Chi magnitude of M = 7.6. Hence, the disturbances in atmosphere and the ionosphere driven by the considered seismic activity have the dimension exceeding 1 Mm. The ionosphere perturbation positioned above the hypocenter of the Chi-Chi EQ in Taiwan will reflect the incident ELF radio waves. We may turn to the F-region ionosphere modifications relevant to the Chi-Chi EQ studied by Chuo et al. (2002) [189] who demonstrated that precursors appear from 1 to 6 days prior to the major shock. We expect that the related modification encompasses the lower ionosphere, while the particular origin mechanism of the ionospheric disturbance is not discussed.
We suggested the following wave interference mechanism for explanation of the Chi-Chi EQ observational data. Two major features of the observed SR modification should be addressed: 1) the single By field component was noticeably increased in the vicinity of the fourth resonance peak and 2) the shift in the peak frequency up to 1 Hz was detected.
We have in mind that global thunderstorm activity was the major source of the ELF field during observations of SR anomaly. The sources are located in the South-East Asia (0˚N and 120˚E), in Central Africa (5˚N and 10˚E), or in the South America (0˚N and 50˚W). We accept a substantial increase in the conductivity of mesosphere and the lower ionosphere which might be treated as a considerable reduction of characteristic heights over the hypocenter of the Taiwan EQ. This kind of ionosphere disturbance will reflect the natural ELF radio waves arriving from thunderstorms. The geometry of the problem is demonstrated in Figure 33 where the propagation paths are shown for the direct and the reflected ELF radio waves. The direct wave arrives from the field source (particular thunderstorm center) and the second wave arises due to wave bouncing from the ionospheric non-uniformity.
Consider initially the Asian thunderstorms (Figure 33). In this case the source-observer distance is equal to 5.5 Mm, while the difference between the geometric path of the direct and the scattered wave is about 5.9 – 5.4 = 0.4 Mm. Obviously, the first constructive wave interference maximum might be expected when this path difference is equal to λ/2 (λ is the wavelength at a given frequency). One readily obtains that λ = 0.8 Mm corresponds to the interference maximum. The SR basic (n = 1) mode wave length is equal to the Earth’s circumference λ = 40 Mm, therefore, the above paths difference might reveal itself at the mode number
, or at the frequency of about 300 Hz. Obviously, the Asian thunderstorms are excluded from the candidates for providing noticeable spectral enhancement around fourth SR mode.
Let us turn to the African and American thunderstorm centers (Figure 33). For Africa, the source-observer distance is equal to 13 Mm. The geometrical path difference is approximately 2 Mm between the direct and the scattered from the EQ waves. For this source, the maximum effect is expected when λ = 4 Mm. The 4 Mm wavelength corresponds to the mode number n = 10, which has the peak frequency of ~60 Hz. Again, we see that the African thunderstorms are more realistic, but not excellent candidates.
The third position of the candidate source is the South America (Figure 33). The source-observer distance in this case is almost antipodal one, about 19 Mm. The geometrical path disparity reaches 4 Mm. The relevant maximum in the constructive wave interference is now shifted at the wavelength λ = 8 Mm. This number corresponds to the expected field maximum at the mode number
![]()
Figure 33. Three thunderstorm regions, (1) Asia (Indonesia), (2) Africa and (3) America (Amazon) and the configurations of the direct path and that scattered at Taiwan are given for each source. After [20] .
or to the frequency of approximately 32 Hz. We may conclude that positioning of the American thunderstorm center at the Amazon basin is the best possible choice to be used in the model computations.
In addition, the geometry of the problem for the Amazon basin is a fortunate one. The great circle arcs from Taiwan and from NAK (Nakatsugawa) toward the Amazon basin are very close, they are almost coincident: the direct (source-observer) and the scattered (EQ-observer) waves have the arrival azimuths separated by ~180˚. This is almost the wave backscatter situation: the direct signal from the American thunderstorms and the signal reflected from the ionosphere non-uniformity over Taiwan will arrive to the observer in Japan along two opposite directions (see Figure 33).
We used three global thunderstorm centers in our computations for verifying that the localized seismogenic ionosphere modification over Taiwan is able to cause an effect similar to experimental observations performed in Japan. Concurrently we evaluate the role of the field source allocation in the emergence of SR anomaly. We do not address here the particular mechanisms responsible for the ionosphere modification. This issue is left for the future work.
The formal approach for solving the ELF scatter from localized ionosphere non-uniformity was described in detail in [190] . The observer is positioned at the NAK observatory in Japan (35.45˚N; 137.3˚E). The knee model [187] is used for the regular conductivity profile in the mesosphere and the lower ionosphere. Its parameters are: the electric characteristic height is close to 55 km and increases with frequency. The magnetic characteristic height is ~95 km and reduces with the frequency increase. These two complex heights allow obtaining the ELF propagation constant ν(f) found from the equation:
, where a is the Earth’s radius and k is the wave number. In the regular cavity, the direct wave from the vertical electric dipole source has the following vertical electric field:
. Here, M(w) denotes the source current moment; Pν(x) is the Legendre function of the complex order ν, and θH is the source-observer angular distance.
The pre-seismic and seismic activity disturbs the conductivity of regular atmosphere in such a way that the undisturbed profile moves downward by 20 km. Formally, this means that the electric height h1 is reduced by 20 km to the disturbed value: hD = h1 – 20 km. The center of ionosphere reduction is found above Taiwan at the point (24˚N; 122˚E). We use in computations the disturbed dimensionless propagation parameter
reaching the maximum value above the center of EQ. The ionosphere perturbation in space
depends on the distance r from its center in the Gaussian manner:
with the horizontal scale factor R = 1000 km.
The electromagnetic problem of wave scattering is solved by using the Stratton-Chu integral equation, which is an equivalent of Maxwell’s equations. The total field in the non-uniform cavity is the sum of the direct E1 and the reflected from the non-uniformity wave E2. The normalized field disturbance is equal to:
, θS is the distance observer-thunderstorm center;
.
Since the size of the localized seismogenic non-uniformity (a few megameters) is much smaller than the wavelength, the ionosphere disturbance might be treated as a Dirac’s delta-function in the surface integral involved in the Stratton-Chu equation.
The physical sense of the above formulas is obvious: the wave reflected from the non-uniformity field is presented as a product of two Green’s functions. The first one accounts for the wave incident at the non-uniformity from the source; this is why its argument is equal to the angular distance (π - θ): source-non-uniformity. The second factor describes the scattered wave itself, it is a function of the observer-non-uniformity distance and the relevant (π - γ) argument. In case of the scalar problem (like reflection of acoustic wave), only the first product remains in the expression for Qν. We treat the vector problem of the electromagnetic fields, and the solution acquires both Green’s functions and their derivatives, such a form accounts for co-existence of electric and magnetic fields.
We demonstrate the results model computations in Figure 34. The frequency variations of the dimensionless field disturbance amplitude
are shown in the left panel. The plots were computed for the point vertical electric dipole sources located in Asia (dotted line), Africa (smooth line), and America (dashed
![]()
Figure 34. Left panel shows the computational results on the frequency dependence for three thunderstorm centres. The right panel indicates the frequency dependences of the vertical electric field expected at our observatory for three sources (from the top to the bottom; America, Africa and Asia). A thin line refers to the direct wave without the effect of ionospheric perturbation, while a thick line, the corresponding result for the disturbed case (with the ionospheric perturbation over Taiwan). After [20] .
line). One may readily observe that the highest normalized disturbance of the SR field is observed for the Asian source. It occurs in the vicinity of 20 Hz frequency or around the third SR mode. It is worth noting that the left frame in Figure 34 demonstrates a “general picture” in the wave scattering and interference. We mean that the high values of the function
might appear owing to a high amplitude of scattered field, however, they may also correspond to rather small amplitudes of the direct wave at an observatory relevant e.g. to the nodal source-observer distance. The lower panel in the right column of plots in Figure 34 indicates that exactly this situation is relevant to the source in Asia.
The three right plots in Figure 34 show the amplitude spectra of vertical electric field component for the Earth-ionosphere cavity excited by the thunderstorms in America, Africa, or on Asia (from top to bottom). The smooth black lines supplied by the “Direct wave” label show the direct ELF radio wave spectra when the ionosphere modification is absent. The dashed lines marked by the “Disturbed” text show the model resonance spectra in these frames when the seismogenic ionosphere perturbation is present. Here, the lower right frame contains the spectra relevant to the Asian thunderstorms. The plots clearly show that the regular field in uniform cavity has a distinct minimum around the 20-Hz frequency. This nodal minimum explains why the large increase was observed in the left frame for the Asian source. The model spectra here show that the Asian source is able to cause an evident disturbance in the vicinity of the second SR mode only.
As it was expected from simplistic geometric consideration, the disturbance driven by the Taiwan seismic activity causes the highest impact on the SR records when the natural sources of ELF electromagnetic radiation are concentrated in South America, specifically, in the Amazon basin. In this case, a two-fold amplitude increase occurs in the vicinity of the fourth SR mode. When the other global thunderstorm centers dominate, a smaller increase is observed in SR amplitude driven by radio wave reflections form the localized non-uniformity positioned over Taiwan. Model computations indicate that the presence of a seismogenic ionosphere perturbation over Taiwan produces the ELF radio wave reflections, which in combination with the direct waves cause a specific interference pattern elevating the amplitude of the fourth SR mode at the NAK observatory in Japan. It should be stressed the concept of wave reflection and subsequent interference effectively explains the major features of experimental SR records: 1) a distinct increase of the fourth SR mode amplitude and 2) the simultaneous displacement of the forth mode peak frequency. As might be observed in Figure 34, the seismo-ionospheric modification may reveal itself at lower frequencies also, but this detail is beyond our goal.
In concluding this sub-section, we must comment on the above used terms ’reflection’ and ’wave interference’. A reader is accustomed to this terminology: it is well familiar and agrees with the intuitive vision of the problem. However, the terms are not exact from the meticulous viewpoint. Strictly speaking, they should be applied when describing an object in the free space, provided that the spatial scales and the characteristic distances noticeably exceed the wavelength. This is not our case. The wavelength at the fourth SR mode is about 10 Mm while the non-uniformity size is about 1 Mm and its distance from the observer is equal to 2 Mm. In addition, we do not deal with the plane incident wave in the free space. Instead, we treat the resonance phenomenon in a closed spherical shell embracing the perfectly conduction globe. Of course, we treated this electromagnetic diffraction problem with maximum precision and performed accurate computations. The results obtained are linked, of course, to the wave diffraction theory. Nevertheless, the results are successfully explained using the conventional physical concept of the wave scattering and interference.
7.4. Recent Advances
In previous sections we addressed the impact of distant (a few Mm) EQs on the SR oscillations. These issues were discussed and generalized recently by Hayakawa et al. (2020a, b, c) [191] [192] [193] [194] . The impact of the nearby (distance less than 1000 km) EQs on SR was separately considered in the paper by Hayakawa et al. (2020c) [193] . The reason was that the clear anomalies in the SR signals were recorded at the NAK field site during two successive and relatively close (~500 km) EQs of the magnitude M » 7. The EQ occurred offshore in the Tohoku province of Japan (the same EQs were already discussed in [124] using the multi-parameter approach). The left side in Figure 35 presents the spectral
![]()
Figure 35. Comparison of experimental (left column) and model (right column) February spectra. After [194] .
anomaly characterized by an increase in the SR amplitude. The modification was detected both prior to and after each of the two nearby EQs. The observed anomaly was interpreted by the seismogenic perturbations of the mesosphere and lower ionosphere. The non-uniformity position was suggested over the EQ hypocenter. Figure 36 illustrates the possible alterations in the vertical profile of atmospheric conductivity addressed in modelling. The regular profile is shown here by the black line. The two disturbed profiles belong either to the compression type with the compression coefficient Kc = 2 (red line) or to the expansion type, Ke = 1.25 (purple line). The details of computations might be found in our papers [191] [192] [193] . The plots indicate that the model results are quite consistent with the observations. An important novelty is connected with the impact of the nearby EQs. The observer-disturbance distance is small in these cases, so that one cannot speak about the wave scattering. The experimental and the model results are explained by the direct positioning of ionosphere depression (elevation) right above the observer. Such a localized height modification alters the excitation factor of the Earth-ionosphere waveguide, which is observed as an increase (or decrease) in the field amplitude. In the theory of radio propagation, the effect is regarded as the impact of the radio wave “take-off-landing platform”.
We have to admit that in the current state of art the most fundamental problem remains unresolved: the absence of reasonable explanation of why and how the ionosphere is modified before and after an EQ. This situation is similar to the condition with the ULF depression, see the discussion in the end of subsection 5.5.
![]()
Figure 36. Vertical profiles of atmosphere conductivity: regular and disturbed (compression type (Kc) and expansion type (Ke)).
7.5. Future Perspectives
The statistical data given in subsection 6.2 illustrate that about thirty observations were linked to the SR anomalies at the NAK field site in Japan when distant (about 2 Mm) EQs occurred in Taiwan. These observations are of a high statistical significance indicating the SR anomalies are the inspiring candidates for the short-term EQ prediction. Therefore, the further case studies are highly desirable combined with the statistical processing of data. The recent cases of detecting the SR anomalies in a close vicinity of EQs (the distance below 1 Mm) indicate the necessity of much more extensive research in this field.
8. Summary and Conclusion
Based on the presentations in this review, we comment on the general remarks on the statistical significance of any physical parameter, which is the main issue of this review. First, we have to say that there are relatively few large magnitude EQ events, which means that most studies of EQ precursors operate on very few examples. An exception is satellite-based studies (e.g., [13] [22] ), but these measurements suffer from resolution issues such as being far from EQ epicenters, and from the measurement locations not being fixed, as already mentioned in this review as well. The lack of large magnitude EQ events also makes it extremely difficult to indicate whether there is a scaling law of effect size with EQ magnitude, although [28] [29] [30] offered some evidence of an effect size more pronounced for larger EQs.
Based on the previous studies over a few decades, it is reached that electromagnetic effects are the most promising candidates as short-term EQ predictors (e.g., [1] - [6] ) and that the ionosphere (both lowest part and also F region) is a plausible main player as short-term EQ prediction. This review has presented a few possible “additional” electromagnetic wave predictors, though not so well investigated in the scientific community.
In addition to the well-established candidate for seismogenic ionospheric perturbations for the short-term EQ prediction, in this review we have suggested additional four-wave phenomena in the ULF/ELF band, including 1) lithospheric ULF radiation (ULF radiation from the lithosphere), 2) ULF horizontal magnetic field depression (as an indicator of lower ionospheric perturbation), 3) ULF/ELF electromagnetic emissions in the atmosphere, 4) SR anomalies (for distant and nearby EQs). As in the previous review [20] , we have presented the essential observational results, with special reference to their statistical significance in terms of probability gain (PB), because we want to estimate the applicability of any physical parameter to real EQ forecast. As a result, we have found that the two wave phenomena of 1) ULF depression effect and 2) ULF atmospheric radiation, are extremely promising with high PG up to 3 or even higher for the real EQ forecast. The next promising is 3) ULF lithospheric radiation, but the number of events is extremely limited, and so the PG by a statistical study for EQs with a magnitude greater than 5 in Kakioka [69] is rather low around 1.5 or so. The most important requirement for the lithospheric ULF radiation will be to accumulate the number of ULF events for large magnitude (M ≥ 7) EQs, together with the application of critical analysis (such as fractal and natural time analyses) (as done in [116] [117] ). The latest phenomenon among the four is (IV) SR anomalies, so the number of events is the poorest as compared with other three-wave phenomena. We are required to accumulate the number of events first of all for huge EQs not only distant but also nearby.
Accumulation of large magnitude EQ events for all four-wave phenomena is the priority, together with the use of critical analysis on the link of the observed anomaly with the relevant EQ. Further, the AI technology will be of great significance in increasing the number of reliable wave events.
With the future directions for the four-wave phenomena, we will be able to obtain more reliable events. Together with this direction, we have to emphasize the importance of multi-parameters (or multi-disciplinary) observations for any big (M ≥ 7) EQs. Even if this kind of multi-parameter is recognized as being of potential use in elucidating the LAIC process [4] [6] , it is a difficult task to accomplish, and so there have been only very few examples in this direction [124] [125] [195] [196] [197] [198] [199] .
Many years ago we established a multi-disciplinary ground-based station in Kamchatka, Russia [200] [201] as a collaborative work between Japan and Russia, being composed of different kinds of observations such as seismic waves, meteorological parameters, geoelectric fields (telluric currents), ULF geomagnetic fields, subionospheric VLF/LF propagation data, etc. The region of Kamchatka, together with the Kurile and the Japanese islands forms the tectonically active northwest margin of the Pacific Ocean. Recently, Chinese colleagues have just established a noble system for monitoring vibrations and perturbations in the LAI (MVP-LAI) in the countryside of Leshan City of Sichuan Province in China in June 2021 [202] . This station looks much more sophisticated than the previous Kamchatka station, by observing 14 physical parameters. Also, this Sichuan region is known to be an EQ-prone region, just like the previous Kamchatka station. Of course, when you establish this kind of coordinated observation for the study of seismogenic effects and elucidation of LAIC process, it is highly required to observe as many candidates for short-term EQ precursors as possible as proposed in this review. In this sense, we are not satisfied with the observing items in this Chinese station, because ULF/ELF wave observation (with sampling of at least 100 Hz) is not included, probably because of the huge amount of data of those ELF waves.
Statistical significance (not only PG, but also any other metrics) of the four ULF/ELF wave phenomena will be highly required to be studied in the future with a more abundant data set. Also, the wave phenomena suggested in this review are recommended to be included as promising EQ predictors for future prospective EQ prediction studies.
Acknowledgements
The 1-second ULF magnetic field data at Kakioka is available from https://www.kakioka-jma.go.jp/obsdata/metadata/en/orders/new/2960, and the source of Kp and Dst index data is the World Data Center for Geomagnetism, Kyoto (https://wdc.kugi.kyoto-u.ac.jp/wdc/Sec3.html).