On the Tempo-Spatial Evolution of the Lower Ionospheric Perturbation for the 2016 Kumamoto Earthquakes from Comparisons of VLF Propagation Data Observed at Multiple Stations with Wave-Hop Theoretical Computations

There have been published many papers on VLF (very low frequency) characteristics to study seismo-ionospheric perturbations. Usually VLF records (amplitude and/or phase) are used to investigate mainly the temporal evolution of VLF propagation anomalies with special attention to one particular propagation path. The most important advantage of this paper is the simultaneous use of several propagation paths. A succession of earthquakes (EQs) happened in the Kumamoto area in Kyusyu Island; two strong foreshocks with magnitude of 6.5 and 6.4 on 14 April (UT) and the main shock with magnitude 7.3 on 15 April (UT). Because the EQ epicenters are not far from the VLF transmitter (with the call sign of JJI in Miyazaki prefecture), we can utilize simultaneously 8 observing stations of our network all over Japan. Together with the use of theoretical computations based on wave-hop theory, we try to trace both the temporal and spatial evolutions of the ionospheric perturbation associated with this succession of EQs. It is found that the ionospheric perturbation begins to appear about two weeks before the EQs, and this perturbation becomes most developed 5 3 days before the main shock. When the perturbation is most disturbed, the maximum change in vertical direction is depletion in the VLF effective ionospheric height of the order of 10 km, and its horizontal scale (or its radius) is about 1000 km. These spatio-temporal changes of the seismo-ionospheric perturbation will be investigated in details in the discusHow to cite this paper: Asano, T. and Hayakawa, M. (2018) On the Tempo-Spatial Evolution of the Lower Ionospheric Perturbation for the 2016 Kumamoto Earthquakes from Comparisons of VLF Propagation Data Observed at Multiple Stations with Wave-Hop Theoretical Computations. Open Journal of Earthquake Research, 7, 161-185. https://doi.org/10.4236/ojer.2018.73010 Received: May 15, 2018 Accepted: July 20, 2018 Published: July 23, 2018 Copyright © 2018 by authors and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/


Introduction
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 [1] [2] [3] [4].
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 [5] [6] [7], but also in the upper F region of the ionosphere [8] [9].Furthermore, such examples are ULF (ultra low frequency)/ELF (extremely low frequency) radiation and ULF depression [10].
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 [3], because such VLF propagation anomalies tend to appear about one week before an EQ.Since the pioneering works by Russian and Japanese [11] [12] [13], 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 [14]), 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 [15] [16] [17] being stimulated by our Japanese network [3] [5].
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 [14].So, statistical studies on the correlation of those VLF/LF anomalies with EQs with M greater than 5.5 or so have been done [7] [18] [19] [20] [21].The latest result by Hayakawa et al.
(2010) [7] 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) [22].
Previous VLF case studies were based on the VLF/LF records for a few relevant propagation paths at maximum.Yamauchi et al. (2007) [26] 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) [27] [28] 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 [31].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.

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 (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.06UT, 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 [31], 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.

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 [3] [5] [6]).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 [12] [13] [21], and the other is "nighttime" fluctuation method [6] [7] [18] [19] [20] [22].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.

Analysis Results
Hayakawa and Asano (2016) [32] 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 [7] [18] [19] [20], 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.(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.

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 [33].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 [34].Hayakawa et al. (1996)  [12] and Molchanov et al. (1998) [13] have used the full wave theory to explain a change in terminator times for the 1995 Kobe EQ, and Molchanov and Hayakawa (1998) [35] have made a statistical study on the correlation between the terminator times and EQs.Later Yoshida et al. (2008) [36] 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. [13].
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 [33].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 R i .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 F i (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 [33] [36]).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) [37] provides the fundamental method, but Wakai et al. [33] 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. [38] 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 [39]).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 h 1 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 h 2 (TX) and h 2 (RX), respectively).In the wave-hop program we can change independently both h 1 and h 2 (TX).Because the ionospheric perturbation appears only close to the TX in our present case, h 2 (TX) is the only parameter for us to alter, while h 2 (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 h 1 is changed (either lower or higher) from the normal value of 90 km at night, and also h 2 (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.

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 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 h 2 (TX) is decreased (going to the right on the abscissa means a decrease in h 2 (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 h 2 (TX)) during subsequent several days.Whereas, the station of AKT exhibits the completely different behavior.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) (h 1 and h 2 (TX)).Of course, it is quite reasonable to assume a larger change in h 1 because of its closeness to the EQ epicenter and a smaller change in h 2 (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 h 1 and h 2 (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 h 2 (TX).The parameter h 1 is fixed to −3 km, −4 km, −5 km, −6 km (correspondingly to the lowering of h 1 by 3, 4, 5 and 6 km) and so on, and we theoretically observe clear depletions for lowering of h 1 by 9 -6 km, which seems to be consistent with the observation on 7 April in the right bottom panel at KTU of    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.

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.(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.Open Journal of Earthquake Research 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  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.

Discussion
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 Open Journal of Earthquake Research 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 [29] [30].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 [5] [7] [18] [19] [20] 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) [32] 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  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)) [40].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 [31].Also, these Kumamoto EQs are the same as the-Kobe 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. [12] and Molchanov et al. [13] 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. [12] 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) [41]: 1) Electrostatic hypothesis based on positive hole carriers (Freund [42]), 2) Chemical channel based on radon emanation (and associated electric field) (Pulinets [43]), 3) Atmospheric oscillation hypothesis (Molchanov and Hayakawa [1], Hayakawa et al. [44]), and 4) Electromagnetic channel (Hayakawa et al. [41]), and Liperovsky et al. (2008) [45] 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 [44].Molchanov et al. (2001) [46] 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) [47], 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 hun-dred 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) [48] 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.

Figure 2 .
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 4 .
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.

Figure 5 .
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.

Figure 6 .
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 h 2 (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 (
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 h 1 .The abscissa means the lowering of h 1 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 h 1 for each path is lowered.Only the propagation paths of JJI-ANA and JJI-IMZ are found to have drastic variations in amplitude when h 1 changes.That is, for smaller changes in h 1 , the amplitude decreases, shows a minimum in amplitude when h 1 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 .
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 h 1 for different propagation paths.Only the h 1 dependences of VLF amplitudes at ANA and IMZ are consistent with the blue curves in the bottom panels.

Fig- ure 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 h 2 (TX) by ≥ 8 km and simultaneous depletion of h 1 by 3 ~ 6 km.

Figure 8 .
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 h 2 (TX) as h 1 as a parameter.

Figure 9
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 h 2 (TX), in which h 1 is decreasing.This figure suggests that the amplitude at the RX does not change a lot with the lowering of h 1 .So that, when we assume that h 1 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, Figure10illustrates 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 h 2 (TX) with a parameter of h 1 (decreasing from the normal height of 90 km with a step of 1 km, a blue arrow indicates the direction of decrease in h 1 ).Judging from this plot, when the height of h 2 (TX) is decreased by 5 km or larger (indicated by a red circle in Figure10), 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 .
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 h 2 (TX) as a function of h 1 .

Figures 11 -
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 .
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.

Figure 12 .
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.

Figure 14 .
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.