Natural Hazards and Earth System Sciences Further study on the role of atmospheric gravity waves on the seismo-ionospheric perturbations as detected by subionospheric VLF / LF propagation

As the target earthquake we have taken a huge earthquake (EQ) named Miyagi-oki earthquake on 16 August 2005 (with magnitude of 7.2) and we have analyzed the 4 month period including the date of this EQ. In addition to our previous analysis on the nighttime average amplitude (trend) and nighttime fluctuation, we have proposed the use of fluctuation power spectra in the frequency rage of atmospheric gravity waves (period=10 min to 100 min) as a third parameter of subionospheric VLF/LF propagation characteristics. Then it is found that this third parameter would be of additional importance in confirming the presence of seismoionospheric perturbations. Finally, we have discovered an important role of lunar tidal effect in the VLF/LF data, which appears one and two months before this large EQ.


Introduction
There have been accumulated a lot of evidences on the presence of precursory seismogenic phenomena not only in the lithosphere, but also in the atmosphere and ionosphere (Molchanov and Hayakawa, 2008).Gokhberg et al. (1989) found the nighttime perturbations in the amplitude (and/or phase) of subionospheric VLF propagation in association with earthquakes (EQs).Hayakawa et al. (1996) then obtained a very convincing evidence of ionospheric perturbations for the Kobe EQ by means of shifts in terminator times Correspondence to: M. Hayakawa (hayakawa@whistler.ee.uec.ac.jp) in the VLF data.So, it is elucidated with the use of subionospheric VLF/LF propagation that the lower ionosphere (D/E regions) is extremely sensitive to an EQ (Hayakawa, 2007).There have been recently published a lot of subionospheric VLF/LF works on this subject by means of (i) case studies (studies for any specific and huge EQs including relatively recent EQs such as Niigata-chuetsu EQ (Hayakawa et al., 2006), Sumatra EQ (Horie et al., 2007a) and so on) and (ii) statistical studies on the correlation between VLF/LF propagation anomalies and EQs (Rozhnoi et al., 2004;Maekawa et al., 2006).By making full use of the advantage of VLF/LF method as integrated measurement, the latter papers have confirmed statistically that there exists a significant correlation between the VLF/LF propagation anomalies (that is, perturbation in the lower ionosphere) and EQs with magnitude greater than 6.0 and with shallow depth (less than 100 km).The lead time is a few days to one week.Later numerous works have followed on the EQ effects onto the upper ionosphere (like F region) by means of bottomside sounding, TEC (total electron contents) measurement etc. (Pulinets and Boyarchut, 2004) and a recent paper by Liu et al. (2007) has concluded that the F region perturbation is statistically correlated with larger EQs, and the lead time is 3-5 days.
Because we understand that the whole ionosphere (not only lower part, but also upper part) is disturbed before an EQ, the next problem we have to solve, is to elucidate the mechanism how and why the ionosphere is disturbed due to the precursory process of the lithosphere.We have already proposed a few possible mechanisms of lithosphereatmosphere-ionosphere coupling (Hayakawa, 2004(Hayakawa, , 2007;;Hayakawa et al., 2006); (i) chemical channel, (ii) acoustic and gravity wave channel and (iii) electromagnetic channel.Simple explanation for each channel is given here.As for the first chemical channel, we expect the emanation of radon or so, leading to a change in atmospheric electricity and then to the re-distribution of ionospheric plasma.As for the 2nd channel, the atmospheric waves play an important role in the coupling, while the electromagnetic wave (like seismogenic ULF and VLF/ELF wave) is the agent to induce an ionospheric perturbation due to the direct heating/ionizaion or due to the particle precipitation from the wave-particle interaction in the magnetosphere in the case of 3rd electromagnetic channel.Molchnov et al. (1993) found that the seismogenic ULF/ELF wave is too weak to induce the ionospheric perturbation in the 3rd channel, so that the remaining two mechanisms (chemical and acoustic and gravity wave channels) should be studied extensively.
There have been reported rather a large number of observational facts in support of the 2nd channel, mainly based on subionospheric VLF/LF propagation data.Molchanov and Hayakawa (1998) studied the terminator times in VLF signal, and found that the modulation of 2, 5, and 10 days (the frequency range of atmospheric planetary waves) is greatly enhanced before an EQ.This finding has led us to conclude that atmospheric oscillations might be involved in this lithosphere-ionosphere coupling.Terminator time (either morning or evening) (Hayakawa et al., 1996) appears once a day, so that the harmonic analysis of terminator times enables us to study the modulation with the order of days.Then, higher frequency modulations have been studied on the basis of nighttime data in subionospheric VLF/LF data.Miyaki et al. (2002) found, for the first time, the enhancements in the VLF fluctuation spectrum in the frequency range of atmospheric gravity waves (AGWs) (the period of 10 min to 1 h), which are well correlated with EQs.Later, Rozhnoi et al. (2004), Shvets et al. (2004), Hayakawa et al. (2005), Horie et al. (2007a), andRozhnoi et al. (2007a, b) have provided more evidences on the importance of AGW effects, and recently Horie et al. (2007b) have obtained a very convincing direct evidence of wave-like structures in the case of Sumatra EQ as a strong support of the AGW effects.Further additional supports have been recently publishied by Muto et al. (2009), Korepanov et al. (2009), and Blaunstein and Hayakawa (2009).These would lend us a further support to the 2nd channel as the lithosphere-ionosphere coupling.
So the purpose of this paper is to deal with a big EQ named Miyagi-oki EQ happened on 16 August 2005 (with magnitude of 7.2 and with depth of 36 km) and to investigate further AGW effects in seismo-ionospheric perturbations.We have already studied the ionospheric perturbations extensively for this EQ by means of ground-based VLF observation (Rozhnoi et al., 2007a;Muto et al., 2009) and also by Demeter satellite observation (Rozhnoi et al.,2007a, Muto et al., 2007).The main emphasis of this paper is to investigate the AGW modulation in VLF/LF data for this Miyagioki EQ and to obtain further evidence on the importance of this AGW influence on seismo-ionospheric perturbations.Finally, an additional new phenomenon of lunar tidal effect is recognized for the first time in our VLF/LF data.

Analysis period, EQs and propagation paths used
We have used the time period of 4 months from 1 June 2005 to the end of September, 2005 including the date of the famous and large Miyagi-oki EQ.This Miyagi-oki EQ happened on 16 August 2005, its magnitude was 7.2 and its depth was 36 km, which is large enough for us to expect any seismo-ionospheric effect (Muto et al., 2009).Table 1 summarizes the EQs (with magnitude greater than 5.5, depth smaller than 100 km) during the above period which might be relevant to the propagation anomalies observed at different observing stations.Date, location, depth, and magnitude of those EQs are indicated in Table 1, based on the information by Japan Meteorological Agency.
Figure 1 illustrates an LF transmitter with call sign of JJY (40 kHz) located in Fukushima prefecture, as indicated by a diamond.The data from three observing stations (Kamchatka (abbreviated as KCK), Moshiri (MSR) and Kochi  (KCH) indicated by stars) are analysed.The EQs in Table 1 are plotted in Fig. 1, and the center of each circle indicates the epicenter, its radius, EQ magnitude and EQ depth is indicated in color.As the wave sensitive area for each propagation path, we have used the 5th Fresnel zone for each path.

LF data analysis method
Figure 2 illustrates the raw data on the observation of amplitude (intensity) of subionospheric LF propagation from the JJY transmitter (40 kHz), which is received at abovementioned three stations in Fig. 1.The left refers to JJY-KCH path, the middle, JJY-MSR and the right, JJY-KCK.The date of Miyagi-oki EQ of our concern is 16 August as indicated by EQ in Fig. 2. The date goes upward from the bottom (1 August) to the top (22 August).As is already known from our previous statistical works (e.g.Maekawa et al., 2006), the average nighttime amplitude (trend) decreases before the EQ, and also the nighttime fluctuation (NF) is greatly enhanced before the EQ.This is qualitatively consistent with our previous finding (Maekawa et al., 2006).
Here we define the two physical quantities for our analysis.We use the data during local night between the terminator time evening and terminator time morning in Fig. 2. First , we estimate the value of dA(t)=A(t)−<A(t)> as shown in the top panel of Fig. 3, where A(t) is the nighttime amplitude at a time t on a particular day (red curve) and <A(t)> is the amplitude at the same time averaged over −30∼−1 day of the current day (blue curve) (all days in the period, are used equally for averageing).By using this residue dA(t) (shown by black curve in the bottom of the top panel of Fig. 3), we define the following two quantities.where N means the nighttime period (UT=12∼16 h (indicated by two grey lines in the bottom of the top panel of Fig. 3) because this period is definitely local night for all of the three propagation paths) In addition to these physical quantities as used before in our many papers (Rozhnoi et al., 2004;Maekawa et al., 2006), we propose an additional use of the AGW influence as a third physical quantity.Figure 3 shows an example of how to estimate the fluctuation spectra.First, we perform the FFT analysis on dA(t) and we obtain the fluctuation spectrum in a wide frequency range of AGW and acoustic wave (AW) in the bottom panel of Fig. 3.We then define dS(f )=S(f )−<S(f )>, in which S(f ) is the fluctuation spectrum for one particular day (red curve) and <S(f )> is the average spectrum averaged over −30 to −1 day of the current day (30 to 1 day before the current day) (blue curve).dS(f ) are indicated by black bars in the bottom of the bottom panel of Fig. 3, and we take only positive dS(f )(dS(f )>0) for our interest and the AGW modulation (M) is defined as follows.

AGW M=
∫ AGW dS (f ) df Range (in frequency) of AGW where dS(f ) is plotted against frequency in Fig. 3, and AGW range is limited by two vertical red lines (left, AGW period (T )=100 m and right, 10 m).
Finally we impose the normalization for these three physical quantities (trend, NF and AGW M), in order to delete any long-term (e.g.seasonal) effects in the following way.current day (one day to 30 days before the current day), and σ is the corresponding standard deviation.First of all we look at the time period just around our target EQ, Miyagi-oki EQ.We look at Fig. 4 for the propagation path of JJY-KCK.Just around the Miyagi-oki EQ, we have found two vertical red bars surrounded by a box in red.A few days before this EQ, we have found that the trend decreases below −2σ line and also the NF exhibits an increase above 2σ line.Together with these, it is found that the 3rd quantity of AGW M index introduced in this paper would offer an additional confirmation on the importance of AGW modulation.Another red bar just after this EQ is also characterized by both a significant decrease in trend (exceeding −2σ ) and increases in NF (exceeding 2σ ) and in AGW M index (well exceeding its 2σ line).Red means that three parameters satisfy the criteria.
Then as seen from Fig. 5 for the next propagation path of JJY-MSR, the trend is found to show a decrease exceeding −2σ and simultaneously the NF exhibits an increase NF exceeding +2σ , both a few days before the EQ.The third parameter, AGW M index is also found to show an increase significantly (well above the 2σ criterion) on the same day.So that, this anomaly one day before the EQ is indicated by a vertical red bar.A few more days before the EQ we find another anomalous peak (very close to the 2σ criteria for both trend and NF, but AGW M index is approaching +2σ criteria).These two anomalous behaviors are considered to be a precursor to this EQ.
Figure 6 for the third propagation path of JJY-KCH is seen to exhibit the similar precursory behaviors as in the former two paths.You can notice two red bars; one is a few days before and another, about one week before the EQ.Red means that all of three physical parameters (trend, NF and AGWM) are found to exceed the corresponding 2σ criteria.Together with these red vertical bars, we notice two yellow bars on the same day of EQ and about one week before the EQ.Yellow means that either two of the three parameters satisfy the threshold.As the conclusion, the VLF propagation is considerably disturbed before the EQ.
Here we describe all of the anomalies seen in Figs. 4, 5 and 6 one by one.First we discuss Fig. 4 for the path of JJY-KCK.We have a look at the 1st box in a dotted blue line and with a yellow vertical bar.The 1st and 2nd parameters (trend and NF) are found to satisfy the corresponding criteria, but the AWG M (as the 3rd parameter) is significantly enhanced, but below 2σ line.Then, this anomaly is likely to be seismogenic, but it is not clear whether this anomaly is related with the small EQ on the same day.The second anomaly is seen in the middle of July.This anomaly is characterized by a decrease in trend, and increases in NF and AGW M. Except the trend, two other parameters are seen to satisfy the 2σ criteria.However, this general tendency is likely to suggest that this anomaly is also seismogenic.The anomaly in the end of July (yellow vertical bar and a broken box in blue) has the same properties (decrease in trend and increase in NF) as the seismogenic effect, but its AGW M index exhibits no enhancement.This means that this peak is not likely to be a seismogenic effect, but due to any other effect.
Next we move on to Fig. 5 for the path of JJY-MSR and Fig. 6 for JJY-KCH.As concerns with the EQ on 19 June we can find red bars on the same day for JJY-MSR and one day after the EQ for JJY-KCH.These might be related with this EQ.Another anomalous peak (yellow bars surrounded by a broken box in blue) for JJY-MSR takes place about one week before this EQ, which may be a precursor to this EQ when thinking of the lead time.Next we have to discuss two broken boxes in blue in the later half of July.The first appearing one seems to be common to the two paths of JJY-MSR and JJY-KCH, while the latter one is seen only for the path of JJY-MSR.When we look at the behaviors of these two boxes, we find that the 1st one seems to be seismogenic because we notice a decrease in trend, an increase in NF and also an increase in AGW M index.On the other hand, the two parameters (trend and NF) follow the same behaviors as for the 1st peak, but we find no response in AGW M index.This means that this peak is not seismogenic, but due to other  4, but some evidence on the lunar effect (with period of 29.5 days).Two anomalous peaks about one month and two months before the main event, as indicated by dotted blue boxes in Fig. 4, are likely to be related with the lunar tidal effect.F-3 means full moon-3 days and so on.
reason, which is already mentioned before in Fig. 4. Now we move on to Figs. 7, 8 and 9, in which we would like to introduce a new idea of the lunar tidal effect in our VLF/LF data.In the field of EQ occurrence studies, few workers have recently suggested an important role of lunar tidal triggering in EQ occurrences (Tanaka et al., 2004;Sue, 2009).This tidal (lunar) effect is found to be very pronounced only for large EQs.First, Fig. 7 indicates that there are present very conspicuous anomalies about one month and about two months before the main Miyagi-oki EQ and these are definitely seismogenic because we find a decrease in trend, an increase in NF and an enhancement in AGW M index.It is interesting to note that these anomaly peaks are incidentally spaced with a period of 28 or 29 days.We know that the lunar period is 29.5 days.The peaks about one month and about two months before the EQ are found to be corresponding to the lunar phase of F-2 days (F means full moon).So, these peaks are considered to be highly due to the lunar effect so that they are the precursory effect of the Miyagi-oki EQ.
Similar peaks are noticed for other propagation paths (such as JJY-MSR in Fig. 8 and JJY-KCK in Fig. 9) for the same Miyagi-oki EQ.Similar anomalies are again to be located in time exactly following the lunar cycle of 29.5 days.The peak about 2 months before the EQ corresponds to F-3 days, while F-3 days is the corresponding lunar phase of the peak one month before the EQ.At any rate, we can infer these anomalies as the pre-seismic conspicuous lunar tidal effect.Tanaka et al. (2004) have found that these lunar effects can be recognized only for larger EQs.This is the reason why we have observed this kind of very regular lunar effect even  in the seismo-ionospheric perturbations because this Miyagioki EQ is large enough to have such lunar tidal effects.

Summary and conclusion
In this paper we have proposed the use of a new and 3rd parameter in analyzing the VLF nighttime data.That is, in addition to the previous two parameters (trend and nighttime fluctuation) the 3rd parameter of AGW modulation (M) index is introduced as defined by the enhancement of AGW fluctuation in the nighttime VLF data.We have already introduced this parameter in our routine-based daily analysis of VLF data.
We have extensively analyzed the AGW M index for a huge EQ of Miyagi-oki EQ, and it is found that this AGW M index is significantly increased exceeding the 2σ criteria for all of the three propagation paths, JJY-KCK, JJY-MSR and JJY-KCH.This fact means that the AGW fluctuations are extremely enhanced for a large EQ, which would lend us a further support to the mechanism, 2nd channel (acoustic and gravity wave channel) in the lithosphere-ionosphere coupling mechanism.Furthermore, the simultaneous increase in AGW M index can be used as a definite confirmation of the presence of seismo-ionospheric perturbations.
We know that there does not happen any EQ even though we have observed significant anomalies (decrease in trend and also increase in NF).When we happen to encounter this kind of anomaly, we have to check whether this is related to any EQ or any other effects (such as geomagnetic storms).During this process, we can say that the use of AGW M index is also important to distinguish between seismogenic and other (such as geomagnetic storm) effects.In this sense, we could find out two anomalous peaks in the nighttime data before the Miyagi-oki EQ; one month and two months before the EQ.The time difference is always ∼29 days, which is the lunar period.For these anomalies, the 3rd parameter, AGW M index is found to be notably enhanced, which suggests strongly that these anomalies are definitely seismogenic.So that, rather than correlating an anomaly with any nearby EQ as described in the previous section, it is much more reasonable and acceptable to attribute these anomalous peaks to the lunar tidal effect.Tanaka and Matsumura (2006) have confirmed the presence of lunar tidal effects in EQ occurrences for this Miyagi-oki EQ.In correspondence with this seismic (mechanical) effect, we have also discovered the first evidence on the same lunar tidal effect in the seismogenic phenomena for this particular EQ.Tanaka et al. (2006), Tanaka and Matsumura (2006) and Sue (2009) have emphasized that this kind of lunar modulation can be seen only for a large EQ.The presence of these lunar effects in our VLF/LF data would suggest the expectation of a huge EQ in one or two months, which would be of valuable information for the short-term EQ prediction.

Fig. 1 .
Fig.1.Locations of our transmitter (wih a call sign of JJY) indicated by a diamond and receiving stations (Moshiri (MSR), Kamchatka (KCK) and Kochi (KCH)) indicated by stars (some other informations are also included, JJI transmitter (diamond) and TYM (Tateyama, Chiba) station, but not used in this study).EQs treated are plotted as circles, for which the center is the epicenter, its radius is proportional to magnitude, and its depth is indicated by color.

Fig. 2 .
Fig. 2. Temporal evolutions in August, 2005 of diurnal variation of signal intensity observed at three stations (left: Kochi (KCH), middle: Moshiri (MSR) and right, Kamchatka (KCK)).Date goes upwards from 1 August (bottom) toward 22 August (top), with the earthquake day indicated by EQ (16 August).Time is given in UT (so that LT in Japan is UT+9 h).

Fig. 3 .
Fig. 3.An example of estimating dS(f ) and the AGW modulation (M) index by means of FFT analysis for dA(t).Top panel indicates the diurnal variation of dA(t) (=A(t)−<A(t)>) (in black) (where A(t) (in red) is the amplitude at a time t on a particular day, and <A(t)> (in blue) is the mean amplitude averaged over −30∼−1 day of the current day) and nighttime part (UT=12-16 h) is used for the estimation of dS(f ).The bottom panel indicates the fluctuation spectrum of nighttime dA(t) (red curve: fluctuation spectrum on the particular day and blue curve, the average spectrum over −30∼−1 day of the current day) and the difference (or residue) dS(f ) is plotted by black bars.The range of AGW is indicated with two limiting values (T =100 m and 10 m as vertical red lines) (together with the range of acoustic waves (AW)), and we take the positive dS(f ) to define the AGW M index.
) 2 dt (only dA (t) ≤ 0) DATA * = DATA−DATA σ where DATA will be one of the three quantities (trend, NF and AGW M), and DATA indicates the corresponding mean value of each quantity averaged over −1∼−30 days of the

Fig. 4 .Fig. 5 .
Fig. 4. Temporal evolutions of three physical parameters (nighttime fluctuation (NF) (top), trend (nighttime average amplitude) (middle), and AGW M index (bottom)) during four months including the date of Miyagi-oki EQ on 16 August 2005 for the path of JJY-KCK.All of the three quantities are normalized by their corresponding standard deviations (σ ) (during −30 to −1 day of the current day).Downward red bars in each panel indicate the times of EQs with length indicating the magnitude.You notice one red box and four blue dotted boxes.In the red box you can find two red vertical bars and red means that all of three physical quantities satisfy their criteria, suggesting that these peaks are seismogenic as a precursor and an after-effect of the EQ.Boxes in dotted blue line mean that there exist an anomaly indicated by vertical yellow bars, but two of the 3 parameters satisfy the 2σ criteria (unsatisfaction of one parameter is shown by a full blue box).Gray zone means the period of no observation.

Fig. 7 .
Fig. 7. Temporal evolutions of the same plot in Fig.4, but some evidence on the lunar effect (with period of 29.5 days).Two anomalous peaks about one month and two months before the main event, as indicated by dotted blue boxes in Fig.4, are likely to be related with the lunar tidal effect.F-3 means full moon-3 days and so on.

Table 1 .
EQs used in this paper, with their characteristics.