There has been an enormous progress on the studies of precursors to earthquakes (EQs) during the last two decades since the 1995 Kobe EQ. It is then found that a lot of evidence of EQ precursors has been accumulated, and especially electromagnetic phenomena are found to be promising candidates for short-term EQ prediction     .
Among different kinds of precursors, there are already a few electromagnetic phenomena, which are found to be statistically correlated with EQs with large magnitude (M) greater than 6 on the basis of long-term data. One typical example is the ionospheric perturbation, which occurs not only in the lower ionosphere    , but also in the upper F region of the ionosphere   . Furthermore, such examples are ULF (ultra low frequency)/ELF (extremely low frequency) radiation and ULF depression  .
The sole method to monitor the former lower ionosphere perturbation is the use of VLF/LF (very low frequency/low frequency) sounding. The amplitude (and/or phase) of subionospheric signals from any VLF/LF transmitters are continuously monitored, and the observed signal parameters are mainly determined by the position of reflection height which depends on the value and gradient of electron density. It is typically 80 km in daytime and is about 90 km at night. This VLF/LF method can provide us with the information on the perturbation in the lower ionosphere, which has been found to be very promising for the short-term EQ prediction  , because such VLF propagation anomalies tend to appear about one week before an EQ. Since the pioneering works by Russian and Japanese    , there have been published many papers on the use of VLF/LF method for the study of the seismo-ionospheric perturbations (see, our latest review  ), and this VLF method is becoming a world trend for short-term EQ prediction as understood from the establishments of VLF networks in different countries    being stimulated by our Japanese network   .
Due to the main advantage of integrated measurement of VLF method in the sense that the observed signal is sensitive to any EQs taking place close to the great circle path between the transmitter (TX) and receiver (RX), it is easy for us to accumulate the number of VLF anomalies  . So, statistical studies on the correlation of those VLF/LF anomalies with EQs with M greater than 5.5 or so have been done      . The latest result by Hayakawa et al. (2010)  has concluded a close correlation between the two on the basis of long-term (7 years) data, so that the presence of precursory VLF anomalies or ionospheric perturbations is evident. A further extensive statistical correlation has also been presented by Rozhnoi et al. (2013)  .
Not only presenting clear statistical evidence on the presence of VLF anomalies, but also the case study on the detailed characteristics of electromagnetic phenomena of each EQ event is important. These case (event) studies can be listed: Tokachi-oki EQ (25 September, 2003, M8.3)   , Niigata-chuetsu EQ (23 October, 2004, M6.8)   , Sumatra EQ in Indonesia (26 December, 2004, M9.0)   , 2011 Tohoku EQ (11 March, 2011, M9.0)   , etc. In these works VLF records are investigated primarily, but we can also look for other related electromagnetic phenomena including atmospheric disturbances and lithospheric effects in order to understand the generation mechanism of seismo-ionospheric perturbations.
Previous VLF case studies were based on the VLF/LF records for a few relevant propagation paths at maximum. Yamauchi et al. (2007)  have compared VLF anomalies for the three propagation paths in Japan, and tried to estimate the spatial extent of the ionospheric anomaly. Also, Horie et al. (2007a, b)   have compared the data at a few Japanese stations from the Australian NWC TX to deduce the wavelike structure of the perturbation in the case of the 2004 Sumatra EQ. The purpose of this paper is that we try to utilize the VLF records for all propagation paths in Japan available at the moment, because the 2016 Kumamoto EQs happened relatively close to the TX with call sign of JJI in Miyazaki prefecture. Eight propagation paths have been extensively utilized to elucidate the spatio-temporal characteristics of the seismo-ionospheric perturbation for the Kumamoto EQs. The detailed study on the spatio-temporal evolution of the seismo-ionospheric perturbation based on comparisons of VLF records at multiple stations with theoretical computations by wave-hop method is the first attempt in these VLF studies. Another valuable point is that the main shock of a succession of the 2016 April Kumamoto EQs is characterized by a huge M = 7.3, which is exactly as big as the disastrous 1995 Kobe  . Finally, we discuss those temporal and spatial variations of the ionospheric perturbation for the 2016 Kumamoto EQs, and compare those with the characteristics for the 1995 Kobe EQ, followed by a brief discussion on the generation mechanism of seismo-ionospheric perturbations.
2. EQs Treated in This Paper
A high seismic activity in Kumamoto was characterized by a succession of three rather huge EQs. The epicenter of the main shock of our interest which happened in Kumamoto (geographic coordinates: 32˚46'55.2''N, 130˚43'33.6''E) on 15 April 2016 (16:25:15.7 UT (universal time); i.e., 16 April 2016 01:25:15.7 JST (Japanese standard time)), with M = 7.3 and a shallow hypocenter (depth ~10 km) as shown as a red dot in Figure 1. It is important to note that about one day before, two strong EQs with M = 6.5 and 6.4 occurred nearly at the same place at a very close epicentral distance, on 14 April 12:26:41.1 UT and 15:03:50.06 UT, respectively. These two EQs are considered to be foreshocks of the main shock. The foreshocks are associated with the Hinaku fault, while the main shock is likely to be related with the nearby Futagawa fault. The M of the main shock was as big as the 1995 Kobe EQ  , and was of the same fault type. These EQs are plotted as a single red dot in Figure 1. So, it seems worthwhile to compare the results for this Kumamoto EQ event with those for the 1995 Kobe EQ.
3. VLF/LF Network in Japanand Analysis Method
We use the following VLF/LF network, which was in operation for the last few years. As for VLF/LF TXs, there are two Japanese TXs: one is located in Ebino, Kyushu with call sign of JJI (frequency = 22.2 kHz) (32.04˚N, 130.81˚E), and the other is JJY (frequency = 40 kHz) located in Fukushima (37.37˚N, 140.85˚E). These two TXs are plotted in Figure 1 with blue dots.
Our network consists of 8 VLF receiving stations all over Japan as indicated in Figure 1 with black dots. From the north they are Nakashibetsu (abbreviated NSB), Suttsu (STU), Akita (AKT), Imizu (IMZ), Katsuura (KTU), Kamakura (KMK), Toyohashi (TYH), and Anan (ANA). All of these stations are equipped with identical receivers that register simultaneously the amplitude and phase of
Figure 1. Relative location of the EQ epicenter (as a red dot), VLF TX (JJI) (as a blue dot) in Kyushu and VLF observing stations (NSB, STU, AKT, IMZ, KTU, KMK, TYH and ANA) (all as black dots). Also, the corresponding great circle path between each VLF observing station and the JJI TX is given by a blue line. And a red circle in Kyushu is the EQ epicenter.
ASK (amplitude shift keying) and MSK (minimum shift keying) narrowband modulated signals in the frequency range of 10 ~ 40 kHz from several TXs. We receive the signals from the above two Japanese TXs (JJY (ASK) and JJI (no modulation)) and three foreign TXs (all American and MSK signals), NWC (Australia), NPM (Hawaii) and NLK (Seattle)) (as shown in    ). The reception is carried out by an electric rod antenna, which measures the vertical electric field component of subionospheric signals. The receiver can record signals with time resolutions ranging from 50 ms to 60 s, but we use a sampling frequency of 1 s for our purpose. Only the amplitude data are only used in the following analysis.
There are two conventional methods of VLF analysis: One is the terminator time method    , and the other is “nighttime” fluctuation method       . In the analysis we will adopt the latter nighttime fluctuation method, in which we use a residual signal of amplitude calculated as a difference between the real signal and the signal averaged during previous one month.
4. Analysis Results
Hayakawa and Asano (2016)  have made the preliminary analysis for this Kumamoto EQ event and the VLF data only for JJI-IMZ path have been studied to suggest the presence of VLF propagation anomaly before the EQ. In this paper, we will perform much more extensive studies on VLF records observed at all stations in Japan, which will enable us to investigate the detailed spatio-temporal evolution of the ionospheric perturbation itself.
The continuous observation with this VLF/LF network has been performed during the last few years, so that we have the data on temporal evolutions of VLF propagation characteristics (nighttime average amplitude) for all propagation paths. Figure 2 illustrates the temporal evolutions of nighttime average amplitude on longer-distance four propagation paths (JJI-STU, JJI-AKT, JJI-KMK and JJI-KTU). Then, the corresponding results for other shorter-distance propagation paths are plotted in Figure 3 (JJI-ANA, JJI-IMZ and JJI-TYH). Unfortunately, the observation at NSB in Hokkaido is missing due to the malfunction of the receiver there just during one week before the EQs, so that we have omitted the observational results at NSB.
The physical parameter plotted in Figure 2 and Figure 3, is the average nighttime amplitude, which is the average of residue amplitude during the nighttime period (JST = 19 h - 29 h). As one datum per day, the nighttime amplitude on a particular day (blue rectangle) is normalized by the standard deviation during the previous 15 days (σ). According to our previous works     , a VLF propagation anomaly is characterized by a decrease in the nighttime VLF amplitude (for example, exceeding −2σ or −3σ) and an enhancement in the fluctuation of VLF signal amplitudes. And, the lead time is, on average, about one week. Only the former parameter of nighttime amplitude is treated in this paper.
Figure 2. Left plot indicates the longer-distance propagation paths (JJI-STU, JJI-AKT, JJI-KMK, JJI-KTU). Right panels illustrate temporal evolutions of the VLF normalized nighttime amplitude for those longer-distance propagation paths during over one month (from March 7 to April 19). The date on the abscissa is given in UT (universal time). The top panel refers to the propagation path from JJI to STU, the second, JJI-AKT, the 3rd JJI-KMK, and the fourth, JJI-KTU. The main shock is indicated by EQ. The horizontal blue line in the data refers to −2σ level and the horizontal red line, to −3σ level. Possible precursors are marked by blue curves with arrows.
Figure 3. Same as Figure 2, but for other shorter-distance propagation paths (JJI-ANA, JJI-IMZ, and JJI-TYH).
Grouping all propagation paths into Figure 2 and Figure 3 has some significance, because Figure 2 refers to the larger propagation distances, while Figure 3 to the smaller propagation distances. A global inspection of Figure 2 and Figure 3 may indicate clear evidence on the presence of VLF propagation anomalies (ionospheric perturbations) in both figures (i.e., depletion of average nighttime amplitude), because the horizontal blue and red lines in the plots refer to −2σ and −3σ respectively. There are seen possible VLF precursors, which are marked by blue curves with arrows. Let us look at each figure carefully, with special attention to the period from about two weeks before the EQ to the EQ time. First of all, we can find some conspicuous depletion in amplitude on the top two panels (JJI-STU and JJI-AKT) in the end of March and in the beginning of April in Figure 2 as indicated by blue curves with arrows. On the other hand, the bottom two panels corresponding to the propagation paths of JJI-KMK and JJI-KTU have exhibited clearer precursory anomalies on those paths with a few anomalies exceeding the −3σ criterion. Namely, on the propagation path of JJI-KMK (third panel on the right), there appears a general tendency of depletion in amplitude (indicated by a blue curve with arrow) starting on 3 April, with the maximum depletion exceeding −3σ on 7 April, and a subsequent period of depletion until 10 April. The date with maximum amplitude depletion is just one week before the 1st foreshock on 14 April. Similarly to the previous panel, the path of JJI-KTU (fourth panel in Figure 2) showed clear depletions during a few days from 29 March to 1 April and on the same day of 7 April as in the 3rd panel (again indicated by a blue curve with arrow).
How about the precursors in Figure 3 with shorter propagation paths? It seems to us that there appear clear tendencies of amplitude depletion for all of these paths (JJI-ANA, JJI-IMZ and JJI-TYH). On the top panel of Figure 3 (JJI-ANA) we have indicated the possible precursory behavior by a blue curve with arrow; that is, the depletion tendency starts on 3 April, with a prolonged period of depletion (exceeding −2σ) during one week till 10 April. Similar tendencies are found for other two paths (JJI-IMZ and JJI-TYH), also indicated by blue curves with arrow. The path of JJI-IMZ exhibited two days of depletion exceeding −2σ, but this tendency is also seen for JJI-TYH (but not so obvious). This kind of similarity is of great importance in finding out VLF precursors, because the simultaneous observation at multiple stations and their comparison as discussed in this paper, is the main advantage of our network observation. However, the amplitude depletion for these short distances seems to have taken place a little earlier than the clear precursor in the previous figure (Figure 2). However, the depletion of 3 April is found to be simultaneous to the one in Fig. 2 on the path to KMK.
5. Some Explanation of Wave-Hop Theory
In order to interpret the phenomena in Figure 2 and Figure 3, we make full use of theoretical considerations. The theoretical method we use in this paper, is so-called “wave-hop” method, which is known to be very effective and useful for relatively short distance (distance less than a few thousand kilometers) propagation at VLF/LF  . The wave-hop method is principally based on ray theory, but it takes into account the wave intensity. So this wave-hop method is situated just between the simple ray theory and the complicated full-wave theory  . Hayakawa et al. (1996)  and Molchanov et al. (1998)  have used the full wave theory to explain a change in terminator times for the 1995 Kobe EQ, and Molchanov and Hayakawa (1998)  have made a statistical study on the correlation between the terminator times and EQs. Later Yoshida et al. (2008)  have utilized the wave-hop theory to account for the terminator time changes, and of course, Yoshida et al’s results are in consistence with the former full-wave result by Molchanov et al.  .
The wave received at a VLF RX is composed of ground wave and sky waves propagated in the Earth-ionosphere waveguide. The computation of sky waves is principally based on the ray theory, and the ground wave is calculated by the full-wave Sommerfeld integrations. More details are given in  . The configuration of our problem is given in the left part of Figure 4, in which you can notice a 1-hop sky wave and a 2-hop sky wave. Generally, we can include many higher-hop waves, such as n-hop wave (n: the number of reflection from the ionosphere). In the case of 1-hop wave, the reflection point (abbreviated as RP in the figures) of the ionosphere is located just in the middle between the TX and RX. The reflection coefficient at the ionosphere is defined by Ri. The ionospheric height is a function of the solar zenith angle (χ), so that the launching elevation angle (ψ) and path length (L) are determined by the distance (d) between the TX and RX and the ionospheric height (h). Though Figure 4 is illustrated for a flat configuration, the real situation is for curved Earth and ionosphere, and these results in an additional effect of focusing expressed by Fi (focusing at the ionosphere). With taking all these effects into account, the electric field at the receiving point can be estimated (see the details in   ).
Figure 4. Leftpanel illustrates the schematic illustration of the use of wave-hop method. TX and RX are the transmitter and a receiver, respectively. One-hop wave reflects from the lower ionosphere at a point just in the middle between the TX and RX (indicated by a red star). On the other hand, the two-hop wave suffers from the ionospheric reflection at two places (indicated by a yellow and a blue star). Right panel refers to an example path from JJI to KUT, with the positions of ionospheric reflections of 1-hop and 2-hop sky waves. RP means reflection point.
Normally it is sufficient to consider only the 1-hop and 2-hop components in our following computations. Recommendation ITU-R (ITU, 2002)  provides the fundamental method, but Wakai et al.  have modified it to a computer-based method. This was followed by a further revision that introduced a reflection height model derived from the parabolic distribution of electron densities in the D and E layers. In this algorithm, different effects of diurnal variation, solar zenith angle, sunspot numbers etc are taken into account. Wakai et al.  have compared the wave-hop prediction with observations in Japan at different distances, and the good agreement between the two suggests the usefulness of this wave-hop technique.
Though the ITU recommendation says that this wave-hop method is preferably used at frequencies above 60 kHz, several papers on the comparison of field intensity by both methods of waveguide and wave-hop, have yielded that a good agreement is obtained between the two, even at frequencies down to 20 kHz (see a recent review by Lynn  ). This is the reason why we adopt a simpler wave-hop method in this paper.
In the right panel of Figure 4 as an example, we consider the propagation from JJI to KTU in Chiba prefecture. The ionosphere just in the middle between TX and RX, is the RP for one-hop sky wave (with h1 as the reflection height). While, there are similarly two reflection points for 2-hop sky waves; that is, one close to the TX, and another, close to the RX (each reflection height defined as h2 (TX) and h2 (RX), respectively). In the wave-hop program we can change independently both h1 and h2 (TX). Because the ionospheric perturbation appears only close to the TX in our present case, h2 (TX) is the only parameter for us to alter, while h2 (RX) is not so important because the ionosphere at the RP of 2-hop wave close to the RX is not expected to be disturbed even if we have a perturbation. Hence, the value of h1 is changed (either lower or higher) from the normal value of 90 km at night, and also h2 (TX) is changed in the following computations.
Figure 5 illustrates the amplitude variation with propagation distance at mid-night (JST = 24 h) at the frequency of JJI transmitter with considering the uniform ionosphere and the resultant amplitude is based on the superposition of one-hop and two hop waves. The TX is JJI at Miyazaki in Kyushu and different observing stations of our network are considered. If any one of the receiving points is located very close to any amplitude sharp minimum due to wave interference, the amplitude there is likely to vary drastically, so that the interpretation of observational data seems to be troublesome. In this sense, Figure 5 suggests that all of the receiving stations are not like this.
6. Theoretical Support by Wave-Hop Theory and Interpretation of Observational VLF Records
We make a reasonable assumption that the ionospheric perturbation begins to appear above the EQ epicenter well before the EQ, and it develops (that is, the degree of perturbation becomes more enhanced together with an enhancement
Figure 5. Theoretical estimation on the VLF amplitude at different distances from the TX(JJI), corresponding to the positions of seven observing stations except NSB. The average nighttime reflection height (90 km) is used, and midnight (JST = 24 h) condition is assumed.
in the spatial extent). Even though we increased the reflection heights, we could not obtain any computational results in consistence with the observation, so that we present the results only when the reflection height goes down. About two weeks before the EQas shown in Figure 6, we have observed, on the bottom right panels, anomalous VLF amplitude changes only on the two propagation paths of JJI-STU and JJI-KTU: amplitude depletion, followed by a gradual recovery (as indicated by a blue curve with arrow). This period seems to be a very initial phase of precursory ionospheric perturbation, which seems to appear very close to the EQ epicenter, probably at the TX side RP of two-hop wave in Figure 4. The solar-terrestrial conditions to be installed in the wave-hop program, were assigned to the time of EQ period. The top right of Figure 6, illustrates the possible amplitude changes expected at all observing stations (with different colors) when h2 (TX) is decreased (going to the right on the abscissa means a decrease in h2 (TX)). EQ in the figure refers to the main shock. No significant changes are expected at ANA, TYH, IMZ and KMK. But you can notice an initial decrease in the theoretical amplitude, and a subsequent increase at the station of KTU. The similar tendency is also theoretically predicted for STU station (also indicated by a blue curve with arrow). This theoretical prediction seems to be consistent with the observational result in the right bottom. This may suggest that the RP of the TX-side of two hop wave begins to decrease and then to develop (that is, more decrease in h2 (TX)) during subsequent several days. Whereas, the station of AKT exhibits the completely different behavior.
Figure 6. Interpretation of VLF data at two stations (STU and KTU in red in the left map) presented on the right bottom panels. Possible precursors are indicated by a blue curve with arrow. The top right panel illustrates the theoretical expectation on the amplitude variation versus h2 (TX) at different observing stations (with different colors), but the two blue curves for KTU and STU in the top right are found to be consistent with the blue curves with arrows in the bottom panels (the explanation of the plots in the right bottom is the same as in Figure 2 and Figure 3, with the blue horizontal line of −2σ and the red horizontal line of −3σ). In the left panel two red stars refer to the RP (STU, 2 hop) and RP (KTU, 2 hop).
Figure 7 (right bottom) illustrates the temporal evolutions observed of VLF amplitude (nighttime average amplitude) for two propagation paths of JJI-ANA and JJI-IMZ, and please look at the period of 4 - 11 April, about 10 days to one week before the EQs. As compared to the previous Figure 6, we can assume that the ionospheric perturbation generated during the period of Figure 6, is going to develop with more enhancement (with more depletion of ionospheric VLF height) and to expand its spatial scale. So that, as given in the left panel we attempted to lower the height of 1 hop wave for the above two propagation paths (JJI-ANA and JJI-IMZ), and the right upper panel illustrates the theoretical amplitudes expected at different stations versus the height of h1. The abscissa means the lowering of h1 up to 10 km from the normal height of 90 km. Different colors in the upper right panel refer to different propagation paths (e.g., TYH, etc) when the height of h1 for each path is lowered. Only the propagation paths of JJI-ANA and JJI-IMZ are found to have drastic variations in amplitude when h1 changes. That is, for smaller changes in h1, the amplitude decreases, shows a minimum in amplitude when h1 is 83 - 82 km (i.e. 7 - 8 km lowering), and then we have found a subsequent increase. The right bottom panels are the observational facts for those two propagation paths, and the anomalies are indicated by
Figure 7. Interpretation of VLF amplitude at two stations of ANA and IMZ presented in the right bottom panels. The tendencies depicted by blue curves with arrow (precursors) are interpreted in terms of the theoretical estimates on the right top panel. The right top panel refers to the variation of VLF amplitude versus h1 for different propagation paths. Only the h1 dependences of VLF amplitudes at ANA and IMZ are consistent with the blue curves in the bottom panels.
blue curves with arrows, which seem to be consistent with the above theoretical expectation.
In the later phase of the previous time period of interest in Figure 7, there is another VLF propagation anomaly only on the propagation path of JJI-KTU as in Figure 8. Please look at the right bottom observational panel which shows one day clear depletion on 7 April on this particular path. During this period, we can expect (or assume) changes of the reflection heights (both for one-hop and two-hop waves) (h1 and h2 (TX)). Of course, it is quite reasonable to assume a larger change in h1 because of its closeness to the EQ epicenter and a smaller change in h2 (TX). Based on this assumption of spatial variation of the ionospheric perturbation, the wave-hop analysis has been performed to estimate the amplitude at KTU, with changing independently h1 and h2 (TX) as in the left panel. The upper right panel of Fig. 8 illustrates the theoretical amplitude expected at the RX (KTU) (on the ordinate) and the abscissa indicates the height of h2 (TX). The parameter h1is fixed to −3 km, −4 km, −5 km, −6 km (correspondingly to the lowering of h1 by 3, 4, 5 and 6 km) and so on, and we theoretically observe clear depletions for lowering of h1 by 9 - 6 km, which seems to be consistent with the observation on 7 April in the right bottom panel at KTU of Figure 8. That is, the observational fact of one day clear depletion on 7 April can be interpreted in terms of a combined effect of depletion of h2 (TX) by ≥ 8 km and simultaneous depletion of h1 by 3 ~ 6 km.
Figure 8. Interpretation of VLF amplitude variation at KTU on 7 April (indicated by a red circle) in the right bottom panel. The top right panel illustrates the wave-hop estimate of VLF amplitude at KTU versus h2 (TX) as h1 as a parameter.
Figure 9 corresponds to our interest in further later time period of 7 April until just before the EQ. The right bottom observational panel illustrates the temporal evolution of VLF anomaly during the above period (encircled by a red circle), which indicates not-so-much change in amplitude for a particular path of JJI-AKT. The right upper panel illustrates the theoretical amplitude expected at the RX (AKT) versus h2 (TX), in which h1 is decreasing. This figure suggests that the amplitude at the RX does not change a lot with the lowering of h1. So that, when we assume that h1 is decreased by more than 5 km, we can expect nearly no change in amplitude, being very consistent with the amplitude change as observed in the right bottom panel.
In the same time period as above, Figure 10 illustrates the situation for a particular propagation path of JJI-STU. The bottom right panel refers to the observational result on VLF amplitude at STU in Hokkaido, and the relevant period of our interest is encircled by a red ellipse. The upper right illustrates the theoretical expectations on the VLF amplitude at STU with the abscissa decreasing the height of h2 (TX) with a parameter of h1 (decreasing from the normal height of 90 km with a step of 1 km, a blue arrow indicates the direction of decrease in h1). Judging from this plot, when the height of h2 (TX) is decreased by 5 km or larger (indicated by a red circle in Figure 10), no significant theoretical change is expected at STU, which seems to be consistent with the observational fact (the nighttime amplitude being just around zero (no significant change)) in the right bottom.
Figure 9. Interpretation of VLF amplitude variation at AKT (encircled by a red circle) with no significant changes in the right bottom panel. The right top panel illustrates the amplitude at AKT versus h2 (TX) as a function of h1.
Figure 10. Same as Figure 9, but for a particular observing station of STU. The area of our interest is indicated by a red ellipse.
7. Summary of Observational Results on Spatio-Temporal Evolution of Seismo-Ionospheric Perturbation
Based on comparisons of the VLF observational results with theoretical works by wave-hop method, we will be able to draw a picture on the temporal and spatial variations of ionospheric perturbations for the 2016 Kumamoto EQs.
Figures 11-16 correspond to the situations of “seismo-ionospheric perturbation” on different days before the main shock. Figure 11 refers to 2 April
Figure 11. The possible area of seismo-ionospheric perturbation associated with the Kumamoto EQ on 2 April (13 days before the main shock (indicated by EQ)) indicated by a vertical line. The right panels are the VLF records at all observing stations. The area with an indication of −3 km means that where the ionosphere is lowered by 3 km.
(13 days before the main shock), and the vertical line indicates the day of our interest, while the right vertical line indicate the date of the main shock (EQ). All of the VLF data at 8 VLF stations are presented on the right of the figure. The large red circle in Kyushu Island is the place of EQs. This Figure 11 indicates that the initial ionospheric perturbation seems to be generated above the EQ epicenter and its perturbation area is indicated by a red area with 3km lowering the ionosphere. Figure 12 corresponds to the day of 4 April (11 days before the main shock), which indicates that the perturbation expands to a larger scale. This figure illustrates an inner circle with ionospheric depletion by 5km, and an outer circle with ionospheric depletion by 3 km. Figure 13 refers to 6 April, 9 days before the main shock, and Figure 14 refers to the day of 8 April (7 days before the main shock). 10 April and 12 April are, respectively, 5 days and 3 days before the main shock. It seems that during this period (5 - 3 days before the main shock) the ionospheric perturbation is likely to be most developed. As in Figure 15 and Figure 16, the ionospheric depletion by 10 km has a radius of 300 km, and that with depletion by 5 km has a radius of 550 km.
Figure 12. Same as Figure 11, but on 4 April (also indicated by a vertical line) (11 days before the main shock). Top possible areas with −5 km and −3 km are given.
One more thing we have to mention, is whether the observed VLF results at some other stations are consistent with the above general tendency of spatio-temporal evolution of the ionospheric perturbation. For example, the anomaly at KMK on 7 and 8 April, VLF change at TYH during 10-12 April, VLF change at ANA during 3 to 10 April, and VLF change at IMZ during 4 to 9 April (as seen in Figures 11-16) are compared with the wave-hop analysis for the above-mentioned spatio-temporal evolution of the perturbation, and we have found that all of these VLF variations observed are consistent with the theoretical speculation based on the above general tendency of ionospheric perturbation.
We can summarize from Figures 11-16 the observational facts together with the help of wave-hop analysis. This paper is the first attempt to investigate both temporal and spatial properties of such a VLF ionospheric perturbation even for a case study, and the following conclusions have emerged from this work.
1) Spatial dynamics
Spatial extension of the ionospheric perturbation in possible association with the 2016 huge Kumamoto EQs can be characterized by two directions (vertical
Figure 13. Same as Figure 11, but on 6 April (also indicated by a vertical line) (9 days before the main shock). An area with −8 km is given.
and horizontal variations). The maximum change in vertical direction is about 10 km: that is, the lower ionosphere responsible for VLF reflection is lowered by about 10 km. This depletion is an effective value, including both the effects of reflection height and density gradient (sharpness). Then the maximum extent of the seismo-ionospheric perturbation may be on the order of 1000 km.
2) Temporal evolution
The seismo-ionospheric perturbation begins to appear about two weeks before the EQ and it is persistent during about 1 week or so. And, the most enhancement in the seismo-ionospheric perturbation appears on 10-12 April, about 5 - 3 days before the main shock, followed by a subsequent decay. As combined with the spatial variations, we can summarize that during the developing phase, the lower ionosphere becomes lower with an average rate of about 1 km per day, while the horizontal extents increases with a rate of 50 - 60 km per day.
The 2016 Kumamoto EQs happened accidentally close to the transmitter of JJI with a distance of 100 km, so that we could utilize the VLF data at all stations of
Figure 14. Same as Figure 11, but on 8 April (also indicated by a vertical line) (7 days before the main shock). The possible areas of −8 km and −5 km are plotted.
our VLF/LF network. The simultaneous use of a multiple of VLF records is the first attempt to estimate the spatio-temporal evolution of the characteristics of lower ionospheric perturbations in possible association with the Kumamoto EQ with the full support of wave-hop theory.
Two foreshocks on the same day happened for the present event, and the main shock took place one day later. Similarly to the case of the 2011 Tohoku EQ, the so-called foreshock with M = 7.3 happened before the main shock with M = 9.0, but we observed the presence of a very clear VLF anomaly on 5 and 6 March, which is considered to be a precursor to the main shock on 11 March   . Hence, in the present Kumamoto event, if any ionospheric perturbation may occur, it is highly likely to be a precursor to the M = 7.3 EQ (main shock) on 15 April (UT).
The previous VLF papers based on the nighttime fluctuation method      have yielded that the average nighttime amplitude has exhibited a significant decrease (below −2σ criterion), together with the enhancement in fluctuation. Also, the lead time is generally about one week. These works have
Figure 15. Same as Figure 11, but on 10 April (also indicated by a vertical line) (5 days before the main shock). The seismo-ionospheric perturbation seems to be most developed.
been and are still effective in order to suggest the presence of precursory ionospheric perturbations, but we cannot obtain any detailed corresponding features (i.e., temporal and spatial) of those ionospheric perturbations, though a larger depletion is found to correspond to a larger M.
Our preliminary study by Hayakawa and Asano (2016)  has indicated only the “presence” of seismo-ionospheric perturbation by using one single propagation path of JJI-IMZ (shown in Figure 3), but no detailed patio-temporal features of the ionospheric perturbation have been performed so far. In the previous papers on VLF anomalies, they have studied only the temporal evolutions of VLF propagation characteristics (amplitude and phase), but those VLF propagation anomalies did not extensively lead to the elucidation of spatio-temporal evolution of the seismo-ionospheric perturbation. So, this is the main purpose of this paper, together with the help of wave-hop method.
The final results of Figures 11-16 illustrate the spatio-temporal evolutions. The VLF records (only the amplitude is used in this paper) at all observing stations (7 stations all over Japan) are simultaneously utilized, together with the
Figure 16. Same as Figure 11, but on 12 April (also indicated by a vertical line) (3 days before the main shock). The ionosphere is still most perturbed on this day, and begins to show a decay afterwards.
help of wave-hop method. It is then found that the seismo-ionospheric perturbation begins to appear above the EQ epicenter about two weeks before the EQ. This perturbation continues to be more developed and most seriously enhanced 5 - 3 days before the main shock, followed by decay. The development of the perturbations means the enhancement in lowering the effective VLF reflection height and also that in the spatial extent. The maximum lowering of the VLF reflection height is about 10 km, and the radius of the perturbation with −5 km reflection height is ~600 km, so that the perturbation seems to extend up to ~1000 km or so. This observational value seems to be consistent with the empirical formula on the EQ preparation radius of R = exp (M) (radius R (in km))  . The rate of lowering the ionosphere is about −1 km/day and that of spatial expansion is ~50 - 60 km/day during the development of the perturbation. This is the first result to obtain the spatio-temporal evolution of the seismo-ionospheric perturbation, as deduced from a combination of VLF records at several stations in Japan and wave-hop method.
This 2016 April Kumamoto EQs had the largest M = 7.3, which is as big as the disastrous 1995 Kobe EQ  . Also, these Kumamoto EQs are the same as theKobe EQ, because both EQ events are of the same fault-type, so that it is worth comparing the results of VLF data (ionospheric perturbation) for the two EQ events. In the case of the 1995 Kobe EQ, its epicenter was located just in the middle of the great-circle path, so that the VLF characteristics obtained are highly likely to reflect clearly the ionospheric perturbation generated. Hayakawa et al.  and Molchanov et al.  have found significant changes in the terminator times, which have appeared from 4 days before the EQ to the EQ date. And, the maximum shifts in terminator times have occurred two days before the EQ. With the help of full-wave method, Hayakawa et al.  deduced the lowering of VLF reflection height by a few kilometers in order to interpret the shift in terminator time. Despite the same M for both Kobe and Kumamoto EQ events, it seems that the perturbation for the Kumamoto EQ event is much more enhanced than the Kobe EQ in the sense of persistence period, enhancement depletion in VLF reflection height, but the general temporal evolutions are much in common such as the maximum perturbation appearing a few days before the EQ. This difference can be explained by the difference in the date of the event. Because the 2016 Kumamoto EQ event happened after the 2011 Mega EQ, the lithosphere all over Japan, even in the western part of Japan such as Kumamoto is likely to be much more disturbed than at the date of Kobe EQ. So it is expected that the corresponding ionospheric disturbance is much more enhanced by an EQ event with the same M after the Mega EQ in 2011 than before it.
Finally, we try to comment on the generation mechanism of seismo-ionospheric perturbations. A few hypotheses have been proposed and summarized in Hayakawa et al. (2004)  : 1) Electrostatic hypothesis based on positive hole carriers(Freund  ), 2) Chemical channel based on radon emanation (and associated electric field) (Pulinets  ), 3) Atmospheric oscillation hypothesis (Molchanov and Hayakawa  , Hayakawa et al.  ), and 4) Electromagnetic channel (Hayakawa et al.  ), and Liperovsky et al. (2008)  made an extensive review on the 2nd ~ 4th channels on the basis of their idea that the characteristics both in space and time before an EQ are very variable which makes us unable to stress only one model and to reject the remaining models. The large spatial scale of the order of 1000 km from the EQ epicenters may be difficult to explain in terms of the 1st and 2nd hypotheses because there should be the radon emanation and crustal stress at such far distances from the EQ epicenter, but it seems to be consistent only with the generation of atmospheric gravity waves (AGWs) propagating very obliquely up to the lower ionosphere  . Molchanov et al. (2001)  provide arguments that water, radon and gas eruptions before an EQ could originate mosaic (in space) and twinkle (in time) spots of atmospheric temperature and density variations leading to the generation of AGW turbulence. As considered by Mareev et al. (2002)  , AGWs produce turbulent variations of density and electric field in the lower ionosphere with a large scale due to oblique propagation effects even though the horizontal scale of the AGW source size on the ground is of the order of ~50―a few hundred kilometers. In the area of Kumamoto where the medium-term EQ prediction indicates that the probability of having M7 class EQs in the coming 30 years is about 1%, there are available only few observations as related to the above lithosphere-atmosphere-ionosphere coupling. Schekotov et al. (2017)  have summarized different electromagnetic effects for this Kumamoto EQ, but no information on radon and crustal changes has been reported up to now. Our future work will be the coordination of different (electromagnetic, mechanical and chemical) data in order to go into details of generation mechanism of the ionospheric perturbation for the 2016 Kumamoto EQ event.
The authors are grateful to Earthquake Analysis Laboratory, which is kind enough to provide us with the VLF data.
 Hayakawa, M. (2011) Probing the Lower Ionospheric Perturbations Associated with Earthquakes by Means of Subionospheric VLF/LF Propagation. Earthquake Science, 24, 609-637.
 Hayakawa, M., Kasahara, Y., Nakamura, T., Muto, F., Horie, T., Maekawa, S., Hobara, Y., Rozhnoi, A.A., Solovieva, M. and Molchanov, O.A. (2010) A Statistical Study on the Correlation between Lower Ionospheric Perturbations as Seen by Subionospheric VLF/LF Propagation and Earthquakes. Journal of Geophysical Research, 115, A09305.
 Liu, J.Y. (2009) Earthquake Precursors Observed in the Ionospheric F-Region. In: Hayakawa, M., Ed., Electromagnetic Phenomena Associated with Earthquakes, Transworld Research Network, Trivandrum, 187-204.
 Liu, J.Y., Chen, C.H. and Tsai, H.F. (2013) A Statistical Study on Seismo-Ionospheric Precursors of the Total Electron Content Associated with 146M ≥ 6.0 Earthquakes in Japan during 1998-2011. In: Hayakawa, M., Ed., Earthquake Prediction Studies: Seismo Electromagnetics, TERRAPUB, 17-29.
 Gokhberg, M.B., Gufeld, I.L., Rozhnoi, A.A., Marenko, V.F., Yampolsky, V.S. and Ponomarev, E.A. (1989) Study of Seismic Influence on the Ionosphere by Super Long Wave Proding of the Earth-Ionosphere Waveguide. Physics of the Earth and Planetary Interiors, 57, 64-67.
 Hayakawa, M., Molchanov, O.A., Ondoh, T. and Kawai, E. (1996) The Precursory Signature Effect of the Kobe Earthquake on VLF Subionospheric Signals. Journal of Communication Research Laboratories, 43, 169-180.
 Molchanov, O.A., Hayakawa, M., Ondoh, T. and Kawai, E. (1998) Precursory Effects in the Subionospheric VLF Signals for the Kobe Earthquake. Physics of the Earth and Planetary Interiors, 105, 239-248.
 Hayakawa, M., Asano, T., Rozhnoi, A. and Solovieva, M. (2018) VLF/LF Sounding of Ionospheric Perturbations and Possible Association with Earthquakes. In: Ouzounov, D., et al., Eds., Pre-Earthquake Processes: A Multi-Disciplinary Approach to Earthquake Prediction Studies, American Geophysical Union, Washington DC, 277-304.
 Biagi, P.F., Piccolo, R., Castellana, I., Maggipinto, T., Ermini, A., Martellucci, S., Bellecci, C., Perna, G.., Capozzi, V., Molchanov, O.A., Hayakawa, M. and Ohta, K. (2004) VLF-LF Radio Signals Collected at Bari (South Italy): A Preliminary Analysis on Signal Anomalies Associated with Earthquakes. Natural Hazards and Earth System Sciences, 4, 685-689.
 Sasmal, S. and Chakrabarti, S.K. (2009) Ionosperic Anomaly Due to Seismic Activities Part 1: Calibration of the VLF Signal of VTX 18.2 kHz Station from Kolkata and Deviation during Seismic Events. Natural Hazards and Earth System Sciences, 9, 1403-1408.
 Raulin, J.P., David, P., Hadano, R.A., Saraiva, C.V., Correia, E. and Kaufman, P. (2009) The South America VLF NETWork (SAVNET). Earth, Moon and Planets, 104, 247-261.
 Rozhnoi, A., Solovieva, M.S., Molchanov, O.A. and Hayakawa, M. (2004) Middle Latitude LF (40 kHz) Phase Variations Associated with Earthquakes for Quiet and Disturbed Geomagnetic Conditions. Physics and Chemistry of the Earth, 29, 589-598.
 Maekawa, S., Horie, T., Yamauchi, T., Sawaya, T., Ishikawa, M., Hayakawa, M. and Sasaki, H. (2006) A Statistical Study on the Effect of Earthquakes on the Ionosphere, Based on the Subionospheric LF Propagation Data in Japan. Annales Geophysicae, 24, 2219-2225.
 Kasahara, Y., Muto, F., Horie, T., Yoshida, M., Hayakawa, M., Ohta, K., Rozhnoi, A., Solovieva, M. and Molchanov, O.A. (2008) On the Statistical Correlation between the Ionospheric Perturbations as Detected by Subionospheric VLF/LF Propagation Anomalies and Earthquakes. Natural Hazards and Earth System Sciences, 8, 653-656.
 Chakrabarti, S.K., Sasmal, S. and Chakrabarti, S. (2010) Ionospheric Anomaly Due to Seismic Activities Part 2: Evidence from D-Layer Preparation and Disappearance Times. Natural Hazards and Earth System Sciences, 10, 1751-1757.
 Rozhnoi, A., Solovieva, M. and Hayakawa, M. (2013) VLF/LF Signals Method for Searching of Electromagnetic Earthquake Precursors. In: Hayakawa, M., Ed., Earthquake Prediction Studies: Seismo Electromagnetics, TERRAPUB, Tokyo, 31-48.
 Shvets, A.V., Hayakawa, M., Molchanov, O.A. and Ando, Y. (2004) A Study of Ionospheric Response to Regional Seismic Activity by VLF Radio Sounding. Physics and Chemistry of the Earth, 29, 627-637.
 Cervone, G., Maekawa, S., Singh, R.P., Hayakawa, M., Kafatos, M. and Shvets, A. (2006) Surface Latent Heat Flux and Nighttime LF Anomalies Prior to the Mω = 8.3 Tokachi-Oki Earthquake. Natural Hazards and Earth System Sciences, 6, 109-114.
 Hayakawa, M., Ohta, K., Maekawa, S., Yamauchi, T., Ida, Y., Gotoh, T., Yonaiguchi, N., Sasaki, H. and Nakamura, T. (2006) Electromagnetic Precursors to the 2004 Mid Niigata Prefecture Earthquake. Physics and Chemistry of the Earth, 31, 356-364.
 Yamauchi, T., Maekawa, S., Horie, T., Hayakawa, M. and Soloviev, O. (2007) Subionospheric VLF/LF Monitoring of Ionospheric Perturbations for the 2004 Mid-Niigata Earthquake and Their Structure and Dynamics. Journal of Atmospheric and Solar-Terrestrial Physics, 69, 793-802.
 Horie, T., Maekawa, S., Yamauchi, T. and Hayakawa, M. (2007) A Possible Effect of Ionospheric Perturbations Associated with the Sumatra Earthquake, as Revealed from Subionospheric Very-Low-Frequency (VLF) Propagation (NWC-Japan). International Journal of Remote Sensing, 28, 3133-3139.
 Horie, T., Yamauchi, T., Yoshida, M. and Hayakawa, M. (2007) The Wave-Like Structures of Ionospheric Perturbation Associated with Sumatra Earthquake of 26 December 2004, as Revealed from VLF Observation in Japan of NWC Signals. Journal of Atmospheric and Solar-Terrestrial Physics, 69, 1021-1028.
 Hayakawa, M., Hobara, Y., Yasuda, Y., Yamaguchi, H., Ohta, K., Izutsu, J. and Nakamura, T. (2012) Possible Precursor to the March 11, 2011, Japan Earthquake: Ionospheric Perturbations as Seen by Subionospheric Very Low Frequency/Low Frequency Propagation. Annals of Geophysics (Italy), 55, 95-99.
 Hayakawa, M., Rozhnoi, A., Solovieva, M., Hobara, Y., Ohta, K., Schekotov, A. and Fedorov, E. (2013) The Lower Ionospheric Perturbation as a Precursor to the 11 March 2011 Japan Earthquake. Geomatics, Natural Hazards and Risk, 4, 275-287.
 Nagao, T., Enomoto, Y., Fujinawa, Y., Hata, M., Hayakawa, M., Huang, Q., Izutsu, J., Kushida, Y., Maeda, K., Oike, K., Uyeda, S. and Yoshino, T. (2002) Electromagnetic Anomalies Associated with 1995 Kobe Earthquake. Journal of Geodynamics, 33, 477-487.
 Hayakawa, M. and Asano, T. (2016) Subionospheric VLF Propagation Anomaly Prior to the Kumamoto Earthquake in April, 2016, Special Papers: The April 2016 M7.0 Kumamoto Earthquake, Japan. New Concepts in Global Tectonics Journal, 4, 273-275.
 Wakai, N., Kurihara, N. and Otsuka, A. (2004) Numerical Method for Calculating LF Sky-Wave, Ground-Wave and Their Resultant Wave Field Strengths. Electronics Letters, 40, 288-290.
 Molchanov, O.A. and Hayakawa, M. (1998) Subionospheric VLF Signal Perturbations Possibly Related to Earthquakes. Journal of Geophysical Research, 103, 17,489-17,504.
 Yoshida, M., Yamauchi, T., Horie, T. and Hayakawa, M. (2008) On the Generation Mechanism of Terminator Times in Subionospheric VLF/LF Propagation and Its Possible Application to Seismogenic Effects. Natural Hazards and Earth System Sciences, 8, 129-134.
 Lynn, K.J.W. (2010) VLF Waveguide Propagation: The Basics. In: Chacrabarti, S., Ed., Propagation Effects of Very Low Frequency Radio Waves, American Institute of Physics, AIP Conference Proceedings 1286, 3-41.
 Hayakawa, M., Molchanov, O.A., NASDA/UEC Team (2004) Summary Report of NASDA’s Earthquake Remote Sensing Frontier Project. Physics and Chemistry of the Earth, 29, 617-625.
 Freund, F. (2009) Stress-Activated Positive Hole Charge Carriers in Rocks and the Generation of Pre-Earthquake Signals. In: Hayakawa, M., Ed., Electromagnetic Phenomena Associated with Earthquakes, Transworld Research Network, Trivandrum, Chapter 3, 41-96.
 Pulinets, S.A. (2009) Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) Model. In: Hayakawa, M., Ed., Electromagnetic Phenomena Associated with Earthquakes, Transworld Research Network, Trivandrum, 235-253.
 Hayakawa, M., Kasahara, Y., Nakamura, T., Hobara, Y., Rozhnoi, A., Solovieva, M., Molchanov, O.A. and Korepanov, V. (2011) Atmospheric Gravity Waves as a Possible Candidate for Seismo-Ionospheric Perturbations. Journal of Atmospheric Electricity, 31, 129-140.
 Liperovsky, V.A., Pokhotelov, O.A., Meister, C.V. and Liperovskaya, E.V. (2008) Physical Models of Coupling in the Lithosphere-Atmosphere-Ionosphere System before Earthquakes. Geomagnetism and Aeronomy, 48, 795-806.
 Molchanov, O.A., Hayakawa, M. and Miyaki, K. (2001) VLF/LF Sounding of the Lower Ionosphere to Study the Role of Atmospheric Oscillations in the Lithosphere-Ionosphere Coupling. Advances in Polar Upper Atmospheric Research, 15, 146-158.
 Mareev, E.A., Iudin, D.I. and Molchanov, O.A. (2002) Mosaic Source of Internal Gravitywaves Associated with Seismic Activity. In: Hayakawa, M. and Molchanov, O., Eds., Seismo-Electromagnetics: Lithosphere-Atmosphere-Ionosphere Coupling, Terra Science Publishing Co., Tokyo, 335-342.
 Schekotov, A., Izutsu, J., Asano, T., Potirakis, S.M. and Hayakawa, M. (2017) Electromagnetic Precursors to the 2016 Kumamoto Earthquakes. Open Journal of Earthquake Research, 6, 168-179.