Driving Source of Change for Ionosphere before Large Earthquake -Vertical Ground Motion-Chia-Hung

: This paper discusses the relationship between the vertical ground motion and ionospheric disturbances before the Kumamoto earthquake on 16 April 2016, in Kyushu, Japan, using the vertical ground motion measured by slant gauges widely distributed in Kyushu, and the NmF2 observed by ionosondes in Japan and another region. We provide evidence that vertical ground motion excites internal gravity waves (IGWs) that disturb changes in the ionospheric plasma density. From the spectral analysis results of the vertical ground motion data, the summation of various period (frequency) components analyzed from the original data of the slant gauge shows a possible correlation with the change of NmF2 before the earthquake. On the other hand, the inﬂuence of the geomagnetic disturbance on vertical ground motion seems to exist. However, we cannot conﬁrm that vertical ground motion is inﬂuenced by the geomagnetic disturbance (Kp index) and that the earthquake is triggered by the geomagnetic disturbance. There are two conditions for the vertical ground motion to disturb variations in the ionospheric plasma density: (1) The effective vertical ground motion period should be shorter than 5 h. In addition, (2) vertical ground motion should continue to exist so that wave energy can be continuously injected into the atmosphere. A possible mechanism with which to modify the ionosphere is discussed. The results of this study can provide a basis for the future ionospheric precursors of earthquakes by using the vertical ground motion.

Ground-based observations help find the time variations but are limited in location.On the other hand, satellite observations depend on their orbit.For example, a satellite in a sun-synchronizing orbit passes over the exact location at the same local mean solar time.For other types of orbits, the region surveyed differs accordingly.Because of these constraints, it is difficult to establish the morphology of the modified ionosphere over the epicenter or even the area of the disturbed ionosphere.Accordingly, the mechanism generating the disturbed ionosphere has yet to be confirmed due to the lack of sufficient observational evidence(s).
A few mechanisms modifying the ionosphere which have been proposed so far are Radon emanation [16], the charge emission from the pressed rocks [17,18], and internal gravity waves [19,20].
The ionospheric disturbances prior to one giant earthquake that occurred in the northern part of Japan at 05:46 UT (14:46 JST) on 11 March 2011 (38 • 6.2 N and 142 • 51.6 E; depth of 29 km; magnitude of M9.0, JMA) have been investigated using the NmF2 observed by ground-based ionosondes and [O + ] measured with the Defense Meteorological Satellite Program (DMSP) F15 satellite at 850 km [21].The ionosonde data were collected from almost all stations worldwide, and the earthquake-related disturbed NmF2 is investigated (the difference between geomagnetic disturbance and earthquake-related NmF2 behavior can be studied using global NmF2).The disturbing driver to modify the ionosphere might be the dynamo electric field (E-field) generated between 100 km and 200 km altitudes, enhancing the eastward E-field during daytime and westward E-field during nighttime.This mechanism explains all ionosphere features that were found during the preparation period of the Tohoku earthquake, such as the equatorward motion of the mid-latitude plasma trough, the enhancement of [O + ] at the equatorward edge of the mid-latitude trough, the enhancement of [O + ] around the geomagnetic equator at the height of DMSP F15, and the behavior of NmF2 in the daytime and nighttime at Khavalovsk, Wakkanai, Kokubunji, and other ionosonde stations.Global ionosphere map (GIM) TEC also shows the TEC increment over the geomagnetic equator during both daytime and nighttime.The idea of an enhanced dynamo E-field also explains the features reported before large earthquakes [9,[12][13][14]22].
The next question is the mechanism of the enhanced dynamo E-field.One possibility is the change of the wind system vertically or horizontally around the dynamo region, which might establish a local dynamo E-field.Another possibility is that microscale neutral wind irregularities reduce conductivity, thus enhancing the E-field produced between the west and east edges of the Earth due to the dynamo action over the epicenter region according to Ohm's law (for ex., [23]).In either case, the main driver might be the internal gravity wave.Although the role of the wave in the F-region ionosphere has been discussed in different kinds of literature [6,6,19,24], we still need direct evidence on the role of the wind.
A further question related to the internal gravity wave is the mechanism and the place of the wave generation.One possible candidate is the small ground motion before a large earthquake, a possibility mentioned by [21].In the past, the possible influence of ground motion on the ionosphere was discussed by [25] based on F-net seismic observations.In this study, we try to find the extremely small ground motion by using NIED Hi-net tiltmeter (slant gauge) data [26][27][28] before the earthquake and the relationship between the ground motions and the ionospheric variations, if it exists.
As the first step to show the role of ground motion on the ionosphere, we discuss in the present paper the foreshock of the 2016 Kumamoto earthquake, which occurred at 21:26 JST (12:26 UT) on 14 April 2016 at 32 • 44.5 N and 130 • 48.5 E with a depth of 11 km and a magnitude of 6.5.After that, the mainshock of the 2016 Kumamoto earthquake occurred at 01:25 JST on 16 April 2016 (16:25 UT on 15 April) at 32 • 45.2 N and 130 • 45.7 E with a depth of 10 km and a magnitude of 7.3.
The space weather parameters which always need to be studied in order to discuss the variation of NmF2 are described in Section 2. The analytical procedure to find the ground motion is described in Section 3. The relation of NmF2 with the Dst, Kp index, and F10.7 is discussed in Section 4.

Space Weather Situation
A study on the effect of large earthquakes on the ionosphere plasma, especially NmF2 and TEC, requires investigations of the following space weather parameters: 3 h geomagnetic disturbance, Kp index, solar radio flux of 10.7 cm radio wavelength, F10.7, and geomagnetic disturbance due to the ring current, Dst index.The dependence of F10.7 on NmF2 was studied by [1].Ref. [29] studied the altitude dependence of the ionospheric electron density on solar activity.The effect of the Dst index on the TEC was systematically studied by [30]. Figure 1 shows these three parameters for two months (1 March-1 May 2016).Three vertical lines mark three earthquakes on 1 April, 14 April, and 16 April.It is noted that F10.7 shows a maximum value during 11-15 April.It might also be noted that the earthquake did not occur when the Kp index showed a maximum of Kp = 6 on 6-7 March during these two months, while the mainshock of the Kumamoto earthquake occurred on 16 April, when Kp was nearly zero.The effect of the magnetic disturbance is discussed later.
geomagnetic disturbance, Kp index, solar radio flux of 10.7 cm radio wavelength, F10.7, and geomagnetic disturbance due to the ring current, Dst index.The dependence of F10.7 on NmF2 was studied by [1].Ref. [29] studied the altitude dependence of the ionospheric electron density on solar activity.The effect of the Dst index on the TEC was systematically studied by [30]. Figure 1 shows these three parameters for two months (1 March-1 May 2016).Three vertical lines mark three earthquakes on 1 April, 14 April, and 16 April.It is noted that F10.7 shows a maximum value during 11-15 April.It might also be noted that the earthquake did not occur when the Kp index showed a maximum of Kp = 6 on 6-7 March during these two months, while the mainshock of the Kumamoto earthquake occurred on 16 April, when Kp was nearly zero.The effect of the magnetic disturbance is discussed later.

Analysis of Hi-Net Data
The NIED Hi-net system provides ground motion data with high sensitivity at ~778 Hi-net stations over Japan.We chose the Hi-net stations along a line by the center of the epicenter of the 2016 Kumamoto earthquake to investigate the ionospheric disturbances caused by the small ground motion before the earthquake.The stations that we investigated are shown in Figure 2.

Analysis of Hi-Net Data
The NIED Hi-net system provides ground motion data with high sensitivity at ~778 Hi-net stations over Japan.We chose the Hi-net stations along a line by the center of the epicenter of the 2016 Kumamoto earthquake to investigate the ionospheric disturbances caused by the small ground motion before the earthquake.The stations that we investigated are shown in Figure 2.
Remote Sens. 2023, 15, x FOR PEER REVIEW 3 of 12 geomagnetic disturbance, Kp index, solar radio flux of 10.7 cm radio wavelength, F10.7, and geomagnetic disturbance due to the ring current, Dst index.The dependence of F10.7 on NmF2 was studied by [1].Ref. [29] studied the altitude dependence of the ionospheric electron density on solar activity.The effect of the Dst index on the TEC was systematically studied by [30]. Figure 1 shows these three parameters for two months (1 March-1 May 2016).Three vertical lines mark three earthquakes on 1 April, 14 April, and 16 April.It is noted that F10.7 shows a maximum value during 11-15 April.It might also be noted that the earthquake did not occur when the Kp index showed a maximum of Kp = 6 on 6-7 March during these two months, while the mainshock of the Kumamoto earthquake occurred on 16 April, when Kp was nearly zero.The effect of the magnetic disturbance is discussed later.

Analysis of Hi-Net Data
The NIED Hi-net system provides ground motion data with high sensitivity at ~778 Hi-net stations over Japan.We chose the Hi-net stations along a line by the center of the epicenter of the 2016 Kumamoto earthquake to investigate the ionospheric disturbances caused by the small ground motion before the earthquake.The stations that we investigated are shown in Figure 2. Figure 3 shows the temporal variation of the vertical component of ground motion measured with a slant gauge at MSIH stations from 1 April to 1 May 2016.MSIH is located 7.79 km east of the epicenter (green marker in Figure 2).The Hi-net observation data sampling rate is 100 points per second.We can obtain the ground displacement per second by integrating the data per second.Then, we resample the data by 1 min to remove the high-frequency components.Figure 3 shows the 30 min running mean variations.The two earthquakes on 14 April and 15 April can be seen as the severe vibration of vertical ground motions at MSIH station.On the other hand, the data further show small disturbances superposed on the baseline of ground motion appear on 1 April, 4 April, 7-8 April, and 11-14 April, which might be precursors of the 14 April earthquake.A pulse was seen on 1 April due to the earthquake which occurred on 1 April (11:39 JST).
and horizontal scales are expressed by geographic longitude and latitudes.The geomagnetic latitude is about minus 10 degrees from the geographic latitude in the Japanese sector.
Figure 3 shows the temporal variation of the vertical component of ground motion measured with a slant gauge at MSIH stations from 1 April to 1 May 2016.MSIH is located 7.79 km east of the epicenter (green marker in Figure 2).The Hi-net observation data sampling rate is 100 points per second.We can obtain the ground displacement per second by integrating the data per second.Then, we resample the data by 1 min to remove the highfrequency components.Figure 3 shows the 30 min running mean variations.The two earthquakes on 14 April and 15 April can be seen as the severe vibration of vertical ground motions at MSIH station.On the other hand, the data further show small disturbances superposed on the baseline of ground motion appear on 1 April, 4 April, 7-8 April, and 11-14 April, which might be precursors of the 14 April earthquake.A pulse was seen on 1 April due to the earthquake which occurred on 1 April (11:39 JST). Figure 4 shows the temporal variations of the power spectral density (PSD) and period (frequency) spectral (top panel), and the total power of the period below 12 h (bottom panel), with 1 min of Hi-net data at the MSIH station.Three earthquakes occurred on 1 April, 14 April, and 16 April, referred to as EQ1, EQ2, and EQ3, respectively.Because of the earthquake ground displacement, this causes the amplitude of the PSD to be very large.In order to highlight the frequency variation of the tiny vertical ground displacement before the earthquake, the upper limit of color in the top panel of Figure 4 is, therefore, set to 0.2, which consequently results in the appearance of saturated values in Figure 4, but this is not the part of this study that is to be investigated.The behaviors of the total power, which are marked as a and b in the bottom panel of Figure 4, are considered to be related to the earthquake (EQ2/EQ3).The common feature is as follows.Before EQ2/EQ3, the PSD appears at the relative maximum on 4 April and 7 April, and then has a large pulse on 14 April (EQ2) and 16 April (EQ3).The PSD variations in the top panel of Figure 4 show that the abovementioned maximum total power before the earthquakes is caused by the amplitude enhancement at the lower period (higher frequency).The main period range on 7 Figure 4 shows the temporal variations of the power spectral density (PSD) and period (frequency) spectral (top panel), and the total power of the period below 12 h (bottom panel), with 1 min of Hi-net data at the MSIH station.Three earthquakes occurred on 1 April, 14 April, and 16 April, referred to as EQ1, EQ2, and EQ3, respectively.Because of the earthquake ground displacement, this causes the amplitude of the PSD to be very large.In order to highlight the frequency variation of the tiny vertical ground displacement before the earthquake, the upper limit of color in the top panel of Figure 4 is, therefore, set to 0.2, which consequently results in the appearance of saturated values in Figure 4, but this is not the part of this study that is to be investigated.The behaviors of the total power, which are marked as a and b in the bottom panel of Figure 4, are considered to be related to the earthquake (EQ2/EQ3).The common feature is as follows.Before EQ2/EQ3, the PSD appears at the relative maximum on 4 April and 7 April, and then has a large pulse on 14 April (EQ2) and 16 April (EQ3).The PSD variations in the top panel of Figure 4 show that the abovementioned maximum total power before the earthquakes is caused by the amplitude enhancement at the lower period (higher frequency).The main period range on 7 April became most evident with periods less than 5 h.This PSD behavior before the earthquake is consistent with previous studies [31][32][33], which reported that the frequencies included in the ground motion shift to a higher frequency toward the earthquake days.
emote Sens. 2023, 15, x FOR PEER REVIEW earthquake is consistent with previous studies [31][32][33], which reported that cies included in the ground motion shift to a higher frequency toward th days.Figure 5 displays the temporal variation in the total power at seven Hi middle of Figure 5 displays the total power of the MSIH station located clos center, while the top and bottom of Figure 5 show the total power of the sta away from the epicenter, as depicted in Figure 1.It is seen that the total p station, except for the HYGH station, increased on 4 April and 7 April prior quakes (EQ2/EQ3).This indicates a trend of the total power towards short frequencies).The spike of total power on 1 April is influenced by EQ1.How power of the HYGH station does not exhibit such a total power spike.Instea total power only displays small periodic oscillations, which is markedly diffe total power of other stations.The reason might be due to the proximity of th tion to the southeast of the epicenter, located near the Pacific Ocean.Afte maximum relative value on 7 April, the total power dropped to a relative mi on 9 April.However, this trend was not noticeable in the YABH, MSIH, an tions, located 27.22 km, 7.28 km, and 32.11 km away from the epicenter, resp total power value did not decrease significantly and remained stable for a the M6.5 earthquake (EQ2) occurred on April 14.We speculate that this feat a precursor to the earthquake.Further data analysis for other earthquakes m conducting because this phenomenon could be one of the features used fo prediction in the future.Figure 5 displays the temporal variation in the total power at seven Hi stations.The middle of Figure 5 displays the total power of the MSIH station located closest to the epicenter, while the top and bottom of Figure 5 show the total power of the stations farthest away from the epicenter, as depicted in Figure 1.It is seen that the total power of each station, except for the HYGH station, increased on 4 April and 7 April prior to the earthquakes (EQ2/EQ3).This indicates a trend of the total power towards short periods (high frequencies).The spike of total power on 1 April is influenced by EQ1.However, the total power of the HYGH station does not exhibit such a total power spike.Instead, the overall total power only displays small periodic oscillations, which is markedly different from the total power of other stations.The reason might be due to the proximity of the HYGH station to the southeast of the epicenter, located near the Pacific Ocean.After reaching its maximum relative value on 7 April, the total power dropped to a relative minimum value on 9 April.However, this trend was not noticeable in the YABH, MSIH, and TMNH stations, located 27.22 km, 7.28 km, and 32.11 km away from the epicenter, respectively.The total power value did not decrease significantly and remained stable for a period before the M6.5 earthquake (EQ2) occurred on April 14.We speculate that this feature might be a precursor to the earthquake.Further data analysis for other earthquakes might be worth conducting because this phenomenon could be one of the features used for earthquake prediction in the future.

Relation among Kp Index, NmF2, and Total Power
We will now discuss the relationship between three parameters, the Kp index, and total power, shown in Figure 6.We are also interested in the possible trigge earthquakes caused by geomagnetic disturbance [34][35][36][37].
The NmF2 observed at Jeju station is shown in Figure 6 as a typical NmF2 be The behavior of the NmF2 is almost the same as that at other ionosonde stations describe in Figures 7 and 8.Although the ionosonde station, Jeju, is located the from the epicenter, it is helpful to study the longitudinal extension of the earthqua noted that the effect of the earthquake extends wider along the longitude than the [9,21].
First, from the trend of the Kp index and the total power in Figure 6a, we can the total power has significantly higher peaks on 19 March, 24 March, 29 March, and 4 April, which indicates that the amplitude of the short period in the spectral o 4 has increased, among which the 1 April one is due to the influence of EQ1.In com with the Kp index, it can be observed that the Kp index was higher during the when the total power was higher, which we call an a-type trend, i.e., a positive r ship.Next, on 7 April, 12 April (two days before the EQ2), and 21 April (a EQ2/EQ3), there were also relatively large peaks of total power.However, the K corresponded to a smaller value, showing an antiphase with the a-type, which wa fied as a b-type change.In summary, the relationship between the Kp index and t power shifts from a-type to b-type during 4 April and 7 April.We speculate that thi shift reflects the ground situation gradually changing toward earthquake day.
The relationship between the NmF2 and the total power shown in Figure 6b similar to that between the Kp index and the total power described above, with two of a-type and b-type.However, NmF2 has some forms that are different from the dex, such as on 19 March, 29 March, and 24 April, which are b-type in NmF2 bu in the Kp index.The relationship between NmF2 and the total power was also fou ing the Tohoku earthquake, which occurred in the north of Japan on 11 March 20 To conclude the relationship description, the total power shows a one-to-one corre ence with NmF2, suggesting that the ground motion influences the ionosphere.

Relation among Kp Index, NmF2, and Total Power
We will now discuss the relationship between three parameters, the Kp index, NmF2, and total power, shown in Figure 6.We are also interested in the possible triggering of earthquakes caused by geomagnetic disturbance [34][35][36][37].The NmF2 observed at Jeju station is shown in Figure 6 as a typical NmF2 behavior.The behavior of the NmF2 is almost the same as that at other ionosonde stations, as we describe in Figures 7 and 8.Although the ionosonde station, Jeju, is located the farthest from the epicenter, it is helpful to study the longitudinal extension of the earthquake.
It is noted that the effect of the earthquake extends wider along the longitude than the latitude [9,21].First, from the trend of the Kp index and the total power in Figure 6a, we can see that the total power has significantly higher peaks on 19 March, 24 March, 29 March, 1 April, and 4 April, which indicates that the amplitude of the short period in the spectral of Figure 4 has increased, among which the 1 April one is due to the influence of EQ1.In comparison with the Kp index, it can be observed that the Kp index was higher during these days when the total power was higher, which we call an a-type trend, i.e., a positive relationship.Next, on 7 April, 12 April (two days before the EQ2), and 21 April (after the EQ2/EQ3), there were also relatively large peaks of total power.However, the Kp index corresponded to a smaller value, showing an antiphase with the a-type, which was classified as a b-type change.In summary, the relationship between the Kp index and the total power shifts from a-type to b-type during 4 April and 7 April.We speculate that this phase shift reflects the ground situation gradually changing toward earthquake day.
The relationship between the NmF2 and the total power shown in Figure 6b is very similar to that between the Kp index and the total power described above, with two trends of a-type and b-type.However, NmF2 has some forms that are different from the Kp index, such as on 19 March, 29 March, and 24 April, which are b-type in NmF2 but a-type in the Kp index.The relationship between NmF2 and the total power was also found during the Tohoku earthquake, which occurred in the north of Japan on 11 March 2011 [38].To conclude the relationship description, the total power shows a one-to-one correspondence with NmF2, suggesting that the ground motion influences the ionosphere.
The readers might be aware that NmF2 steeply decreases right after 16 April.We cannot explain this rapid reduction only as the reduction of F10.7 [36].One reason might be due to the seasonal dependence of NmF2, and one more possible factor is the effect of many small earthquakes which occurred after 16 April.The readers might be aware that NmF2 steeply decreases right after 16 cannot explain this rapid reduction only as the reduction of F10.7 [36].One rea be due to the seasonal dependence of NmF2, and one more possible factor is th many small earthquakes which occurred after 16 April.

Latitude Variation of NmF2
The effect of large earthquakes on ionosphere NmF2 varies depending on sity of the earthquake-related electric field and geomagnetic latitude.Regardin gitudinal extension, the effect of earthquakes on the ionosphere extends to a mu distance [39].For example, the deviation of plasma drift from ambient plasm tends to a total of 120 degrees in longitude and 30 degrees in latitude for the e that occurred on 16 October 1981 with the epicenter located at −33.134°N and with a magnitude of 7.2, and a depth of 33 km.The latitudinal extension coin the value calculated according to the equation proposed by [40]

Latitude Variation of NmF2
The effect of large earthquakes on ionosphere NmF2 varies depending on the intensity of the earthquake-related electric field and geomagnetic latitude.Regarding the longitudinal extension, the effect of earthquakes on the ionosphere extends to a much greater distance [39].For example, the deviation of plasma drift from ambient plasma drift extends to a total of 120 degrees in longitude and 30 degrees in latitude for the earthquake that occurred on 16 October 1981 with the epicenter located at −33.134 • N and −73.074 • E, with a magnitude of 7.2, and a depth of 33 km.The latitudinal extension coincides with the value calculated according to the equation proposed by [40].Accordingly, we studied the NmF2 behavior observed at four ionosonde stations in Japan and one in Korea.These stations, starting from low latitude, are Okinawa (26. Figure 7a-e show the NmF2 at the five stations from 1 March to 31 March to describe the behavior of NmF2 related to the earthquake that occurred on 1 April (EQ1).It is noted that the daytime NmF2 at Okinawa and Wakkanai show a relatively large value on 23 March, while that at Yamagawa, Jeju, and Kokubunji show a relatively small value.The large NmF2 at Okinawa is due to the plasma lifted slightly by the eastward E-field, while the large NmF2 at Wakkanai is due to the downward plasma flow along the magnetic line of force.The plasma is lifted widely over the epicenter region by the daytime eastward enhanced E-field, and the flux tube of the geomagnetic field is partly filled by the plasma, which is lifted.The small NmF2 at Yamagawa, Jeju, and Kokubunji is caused by an uplift of the plasma to a much higher altitude caused by the eastward E-field that is more substantial than that at Okinawa.It is noted that the variation of ionospheric plasma density is affected by geomagnetic conditions.Figure 1 shows that the average Kp index on 23 March is 2, indicating that 23 March is a geomagnetically quiet day.The variations in NmF2 at the abovementioned stations can be seen as a reference for the geomagnetic calm days.
During 26-31 March, the NmF2 at Okinawa, Yamagawa, Jeju, and Kokubunji showed a relatively small value on 29 March, while that at Wakkanai showed a relatively large value.The small NmF2 at the four stations suggests that the eastward E-field became more substantial than that on 23 March.Daytime NmF2 at Wakkanai shows a larger NmF2 than on 23 March.Nighttime NmF2 shows a more significant value and gradually decreases.The feature suggests that the daytime E-field on 29 March is more enhanced than that on 23 March.According to Figure 6a, the Kp index and total power both show an a-type relationship on 23 and 29 March.However, Figure 6b shows that the relationship between the NmF2 and total power is a-type on 23 March, but b-type on 29 March, meaning that the decrease in NmF2 is associated with an increase in total power.Therefore, the enhanced E-field on 29 March might be caused by the power spectral change of ground motion.It might be seen as evidence of the correlation between the power spectral density of vertical ground motion and the maximum density NmF2.
Figure 8 is shown to discuss the effect of earthquakes EQ2 and EQ3.Unlike Figures 7  and 8 showing a relatively small value around 12 April for all daytime NmF2 at Okinawa, Yamagawa, Jeju, and Kokubunji, in comparison, daytime NmF2 at Wakkanai is shown to be slightly increased on 13 April.The nighttime NmF2 shows a value more significant than that for the earthquake on 1 April, just after the dynamo E-field changes the direction from eastward to westward.The feature suggests that the dynamo E-field for EQ2 is more enhanced than that of EQ1.

Discussion
We report the possible correlation between the power spectral density of vertical ground motion and the maximum density of NmF2 to explain the behavior of NmF2 by the internal gravity wave driven by the ground motion.The ground motion that lasts at least one day produces the atmospheric wave, amplified as it propagates to a higher altitude.At the height of around 90 km, the amplitude is intensified up to ~10 6 times because the neutral density decreases to 10 13 particles/cc compared with the ground 10 19 particles/cc.The amplified waves might produce two phenomena.One phenomenon is the formation of a local wind system.Another phenomenon is that the waves reaching around 90 km start breaking and producing neutral density irregularities, finally producing plasma irregularities.The electric conductivity decreases due to the irregularities, so the E-field enhances during the day and night [41].During the daytime, the eastward dynamo E-field enhances.When the E-field is not intensified but is still stronger than that during the non-earthquake preparation period, the NmF2 around the epicenter region increases because of the lower recombination of electrons with ambient neutral particles, while the increase of the NmF2 is challenging to identify.Further, the intensified E-field lifts the ionosphere to higher altitudes, which reduces NmF2, because plasma is redistributed at a wider height region.The above statements explain the increase or reduction of NmF2, seen in the Okinawa, Yamagawa, Jeju, Kokubunji, and Wakkai ionosonde stations.
During the nighttime, plasma lifted during the daytime is pushed down.If the plasma pushed down is insufficient, both NmF2 and h'F decease, usually seen for earthquakes in the equatorial region.At high latitudes where the magnetic field line is connected to the region over the epicenter, the NmF2 during the daytime shows a slight increase because the plasma lifted to a higher altitude over the epicenter flows down along the geomagnetic line of force.During the nighttime, the plasma, which is accommodated in the magnetic flux tube, flows down along the magnetic line of force.This mechanism increases the nighttime NmF2.If the plasma lifted during the daytime is insufficient, the NmF2 increases right after the electric field changes direction from east to west.However, with the plasma, which was stored in the magnetic tube of flux, the NmF2 reduces.This explains the behavior of the NmF2 at Wakkanai.The effect of the E-field on the behavior of NmF2 is roughly discussed by [42].
The two-day oscillation of NmF2 is well-noticed for Wakkanai NmF2 from 10 April before the earthquake.The frequency spectral analysis of the NmF2 at Kokubunji shows the same.This might result from the interaction between planetary waves and earthquakerelated atmospheric waves [36].

Conclusion and Future Tasks
This paper first reports the direct evidence of the change of NmF2 caused by the small ground motion before the 2016 Kumamoto earthquake.As the same feature appeared for the 2011 Tohoku-Oki earthquake on 11 March 2011, we believe the feature found for the Kumamoto earthquake is relevant.Two conditions for generating the continuous abnormal behavior of ionosphere NmF2 are small ground motion continuously injecting the energy to the atmosphere, and the period of the variation of ground motion being less than 5 h.It is also noted that the 4-5 h of daytime variation of NmF2, which is caused by more extended period waves, is also well-recognized.
However, we concluded that an original driver to modify NmF2 is a periodic fluctuation contained in a small ground motion, not the Kp index, because sometimes the increase of Kp does not coincide with the increase of PSD.This issue must be pursued for other earthquakes that did not occur during the magnetic disturbance.In the future, we will add more earthquake events to verify the relationship between vertical ground motion and ionospheric plasma density disturbances before an earthquake.

Figure 1 .
Figure 1.Temporal variations of solar flux index F10.7,geomagnetic disturbance indices of Dst, and Kp index from 1 March to 1 May 2016.Three vertical lines inserted in the three panels are the days of three earthquakes which occurred at 11:39 JST (02:39 UT) on 1 April with a magnitude of 6.5, at 21:26 JST (12:26 UT) on 14 April with a magnitude of 6.5, and at 01:25 JST on 16 April (16:25 UT on 15 April) with a magnitude of 7.3, respectively.

Figure 1 .
Figure 1.Temporal variations of solar flux index F10.7,geomagnetic disturbance indices of Dst, and Kp index from 1 March to 1 May 2016.Three vertical lines inserted in the three panels are the days of three earthquakes which occurred at 11:39 JST (02:39 UT) on 1 April with a magnitude of 6.5, at 21:26 JST (12:26 UT) on 14 April with a magnitude of 6.5, and at 01:25 JST on 16 April (16:25 UT on 15 April) with a magnitude of 7.3, respectively.

Figure 1 .
Figure 1.Temporal variations of solar flux index F10.7,geomagnetic disturbance indices of Dst, and Kp index from 1 March to 1 May 2016.Three vertical lines inserted in the three panels are the days of three earthquakes which occurred at 11:39 JST (02:39 UT) on 1 April with a magnitude of 6.5, at 21:26 JST (12:26 UT) on 14 April with a magnitude of 6.5, and at 01:25 JST on 16 April (16:25 UT on 15 April) with a magnitude of 7.3, respectively.

Figure 2 .
Figure 2. Locations of seven NIED Hi-net (slant gauge) stations and the epicenter of the 2016 Kumamoto earthquake.A yellow triangle symbol with a black edge marks the epicenter.The vertical and horizontal scales are expressed by geographic longitude and latitudes.The geomagnetic latitude is about minus 10 degrees from the geographic latitude in the Japanese sector.

Figure 3 .
Figure 3.The temporal variation of the vertical component was measured by the slant gauge station of MSIH station from 1 April to 1 May 2016.The horizontal axes are scaled by JST (Japan Standard Time (JST) = UT + 9 h) for the slant gauge signal.

Figure 3 .
Figure 3.The temporal variation of the vertical component was measured by the slant gauge station of MSIH station from 1 April to 1 May 2016.The horizontal axes are scaled by JST (Japan Standard Time (JST) = UT + 9 h) for the slant gauge signal.

Figure 4 .
Figure 4.The temporal variations of the power spectrum density (top) and the total p at MSIH station from 23 March to 22 April 2016.There were three earthquakes occ and 16April, which referred to as EQ1, EQ2, and EQ3, respectively.a and b refer to t two highest PSDs before EQ2/EQ3.In order to see the small frequency variation befor the upper limit value of color in top panel is set to 0.2.This is the reason why the sa occurred in the PSD.

Figure 4 .
Figure 4.The temporal variations of the power spectrum density (top) and the total power (bottom) at MSIH station from 23 March to 22 April 2016.There were three earthquakes occurred on 1, 14, and 16 April, which referred to as EQ1, EQ2, and EQ3, respectively.a and b refer to the times of the two highest PSDs before EQ2/EQ3.In order to see the small frequency variation before earthquakes, the upper limit value of color in top panel is set to 0.2.This is the reason why the saturated values occurred in the PSD.

Figure 5 .
Figure 5. Daily temporal behavior of the total power at seven different stations from 23 Ma April.The MSIH station is closest to the epicenter, while the IMRH and the HYGH station from the epicenter (see Figure 1 to know the distance from the epicenter).

Figure 5 .
Figure 5. Daily temporal behavior of the total power at seven different stations from 23 March to 22 April.The MSIH station is closest to the epicenter, while the IMRH and the HYGH stations are far from the epicenter (see Figure 1 to know the distance from the epicenter).

Figure 6 .
Figure 6.Kp index, NmF2 at Jeju, and total power at IMRH from 1 March 2016 to 31 April (scaled by JST).(a) displays the comparison of the Kp index (black bar) and the total power (red line), while (b) displays the comparison of the NmF2 (black line) and the total power (red line).The marks of 'a' and 'b' indicate the different type of relationship which is described in the text.Three earthquakes occurred on 1, 14, and 16 April, which referred to as EQ1, EQ2, and EQ3, respectively.The mark of '#' in y-axis refers to electrons.

Figure 6 .
Figure 6.Kp index, NmF2 at Jeju, and total power at IMRH from 1 March 2016 to 31 April (scaled by JST).(a) displays the comparison of the Kp index (black bar) and the total power (red line), while (b) displays the comparison of the NmF2 (black line) and the total power (red line).The marks of 'a' and 'b' indicate the different type of relationship which is described in the text.Three earthquakes occurred on 1, 14, and 16 April, which referred to as EQ1, EQ2, and EQ3, respectively.The mark of '#' in y-axis refers to electrons.

Figure 6 .
Figure 6.Kp index, NmF2 at Jeju, and total power at IMRH from 1 March 2016 to 31 Apri by JST).(a) displays the comparison of the Kp index (black bar) and the total power (red lin (b) displays the comparison of the NmF2 (black line) and the total power (red line).The ma and 'b' indicate the different type of relationship which is described in the text.Three eart occurred on 1, 14, and 16 April, which referred to as EQ1, EQ2, and EQ3, respectively.The '#' in y-axis refers to electrons.

Figure 7 .
Figure 7. NmF2 (MHz) at five stations from 1 March to 31 March.The black-dashed line is t of the daily maximum NmF2 value.The blue and green box indicate the date of 23 and 2 respectively.The mark of '#' in y-axis refers to electrons.

Figure 7 .
Figure 7. NmF2 (MHz) at five stations from 1 March to 31 March.The black-dashed line is the curve of the daily maximum NmF2 value.The blue and green box indicate the date of 23 and 29 March, respectively.The of '#' in y-axis refers to electrons.

Figure 8 .
Figure 8. NmF2 at five stations during the period of 1 April to April.The black-dashe curve of the daily maximum NmF2 value.The blue box indicates the date of 12 April.T '#' in y-axis refers to electrons.

Figure 8 .
Figure 8. NmF2 at five stations during the period of 1 April to 30 April.The black-dashed line is the curve of the daily maximum NmF2 value.The blue box indicates the date of 12 April.The mark of '#' in y-axis refers to electrons.