Sub-ionospheric VLF signal anomaly due to geomagnetic storms: a statistical study

We investigate quantitatively the effect of geomagnetic storms on the sub-ionospheric VLF/LF (Very Low Frequency/Low Frequency) propagations for different latitudes based on 2-year nighttime data from Japanese VLF/LF observation network. Three statistical parameters such as average signal amplitude, variability of the signal amplitude, and nighttime fluctuation were calculated daily for 2 years for 16–21 independent VLF/LF transmitter–receiver propagation paths consisting of three transmitters and seven receiving stations. These propagation paths are suitable to simultaneously study high-latitude, low-mid-latitude and midlatitude D/E-region ionospheric properties. We found that these three statistical parameters indicate significant anomalies exceeding at least 2 times of their standard deviation from the mean value during the geomagnetic storm time period in the high-latitude paths with an occurrence rate of anomaly between 40 and 50 % presumably due to the auroral energetic electron precipitation. The mid-latitude and lowmid-latitude paths have a smaller influence from the geomagnetic activity because of a lower occurrence rate of anomalies even during the geomagnetically active time period (from 20 to 30 %). The anomalies except geomagnetic storm periods may be caused by atmospheric and/or lithospheric origins. The statistical occurrence rates of ionospheric anomalies for different latitudinal paths during geomagnetic storm and non-storm time periods are basic and important information not only to identify the space weather effects toward the lower ionosphere depending on the latitudes but also to separate various external physical causes of lower ionospheric disturbances.


Introduction
Very Low Frequency (VLF) signals (3-30 kHz) from the powerful transmitters around the world are an important remote sensing tool to monitor near geo-space environment.If VLF signals are received beyond the skip zone of a transmitter, the amplitude and phase characteristics of received waveforms are very sensitive to the lower ionosphere conditions around transmitter-receiver great circle path (GCP) due to changing propagation characteristics in the Earth-Ionosphere WaveGuide (EIWG) where the VLF waves propagate.Various causes of disturbances have been identified mainly from space and tropospheric phenomena.Therefore VLF signals contain information about an evolution of the lower ionosphere (e.g., D-layer formation and disappearance time) (Chakrabarti et al., 2010), space weather effects such as geomagnetic storms (Kikuchi and Evans, 1983), effects of solar flares and solar energetic particles events (Mitra, 1974;Thomson and Clilverd, 2001;Todoroki et al., 2007), lightning-induced energetic particle precipitations and direct heating due to intensive lightning discharges and Transient Luminous Events (TLEs) (Inan et al., 1985(Inan et al., , 1988;;Hobara et al., 2001) and even effects of strong gamma ray bursts far away beyond the solar system (Tanaka et al., 2010).
One of the most recent but highly challenging and important tasks for the humanity is to study the VLF propagation related to the possible effects of electromagnetic coupling between physical processes in the lithosphere and overlaying ionosphere preceding major seismic activity because of the possible earthquake prediction in the future to mitigate casualties (see, e.g., Hayakawa and Hobara, 2010).
The most problematic issue upon using VLF signals to identify pre-seismic ionospheric anomaly is that received VLF signals are influenced by a large number of physical phenomena.This means that there are many anomalies that may not be associated with seismic origin.Therefore as a first step every possible physical cause affecting VLF signal should be identified and quantitatively investigated to separate the disturbances from seismic origin.Then the coupling mechanism between Lithosphere, Atmosphere, and Ionosphere (LAI) will be comprehensively investigated.
The present paper is devoted to the statistical study of VLF anomalous associated with one of the most pronounced agents to disturb the lower ionosphere; the space weather disturbances as geomagnetic storms.This study is based on the data from 27-month measurements of 18 different transmitter-receiver paths consisting of seven different receiving stations in Japan from three overseas VLF transmitters.
As a result high-latitude paths are found to be more sensitive to the ionospheric perturbations from geomagnetic activities rather than those of low latitude.

Data and instrumentation
Amplitude and phase of VLF communication transmitter signals with emitting frequencies around 20 kHz are continuously monitored by seven receiving stations (NSB, MSR, CHF, MYK, KSG, TYM and KCH) in Japan to study the ionospheric perturbations in the D-layer.Geophysical coordinate and transmitting frequency of VLF stations are summarized in Table 1.Geographical distributions of our VLF In this work, we used the VLF data for the time period between January 2012 and March 2014.During the entire time period of analysis, we used six receiving stations except Moshiri (MSR) because of the data availability.For our data analysis purpose, diurnal variations of electric amplitude from every VLF transmitter-receiver pair with a 2-min temporal resolution were utilized.Moreover we chose the We used the hourly values of Dst and AE indices to specify the time periods of space weather effects.Both data are available from the World Data Center (WDC) for geomagnetism, Kyoto.The Dst index is derived from a network observations of the geomagnetic field located near the equatorial region and is used to identify the geomagnetic storm time period in this study.The AE index is derived from a number of stations located in the auroral zone to monitor the auroral electrojet and is used here to determine the geomagnetically quiet time period (no magnetospheric substorms).

Data pre-processing methodology
The ionospheric VLF signatures of nighttime D-region are extracted from every VLF path from its diurnal amplitude variations.We use nighttime VLF amplitude for data analysis because ionospheric perturbations due to space weather effect such as geomagnetic storms are pronouncedly observed in nighttime rather than daytime (Araki, 1974;Kleimenova et al., 2004).Moreover, the observed VLF amplitude is stronger in general for nighttime due to smaller absorptions during the propagation unless amplitude decreases significantly due to the mode interference.In this study nighttime was defined as the time period during the complete darkness for every transmitter-receiver path.
The nighttime VLF amplitude time series were transformed to three daily statistical parameters namely Trend, Dispersion and Nighttime Fluctuation (NF).These parameters have been considered to be useful to identify the days for strong VLF propagation anomaly due to ionospheric perturbations, which is so-called "nighttime fluctuation method" (Rozhnoi et al., 2004).In this method, we estimated the residual of the signal amplitude dA(t) defined as dA(t) = A(t) − A(t) , where A(t) is the signal amplitude at a time t on a particular day.A(t) is the average amplitude at the same time t over the previous 15 days (i.e., ensemble mean) (in this work) from the current day.This process tends to enhance the short-term variations and reduce the long-term variations such as seasonal dependence.
The first parameter, Trend is the change on the daily nighttime average amplitude.The Trend is obtained daily by averaged over the nighttime period for each day dA(t).The second parameter, Dispersion is based on the conventional statistical quantity of standard deviation (square-root of variability) of dA(t) for each day and is also obtained daily.The third parameter, NF is defined by the area of dA(t) < 0 integrated over the whole nighttime.In fact, this method has been applied extensively to identify the ionospheric anomaly related to the seismic activity as significant decrease in Trend associated with increase in both Dispersion and NF (e.g., Rozhnoi et al., 2004;Hayakawa et al., 2011).

VLF anomalies associated with geomagnetic storm
Diurnal variations of VLF amplitude for high-latitude (NSB-NLK), mid-latitude (NSB-NPM) and north-south (NSB-NWC) paths with Dst values for 16 days around the geomagnetic storms (including the dates before and after the storm) are shown in Fig. 2. The regular diurnal patterns of the VLF amplitude are recognized in all three transmitter signals NLK (top panel), NPM (second panel) and NWC (third panel) before 7 March (the day of the geomagnetic storm onset) with a sharp drop of the Dst index (bottom panel).In general, amplitude of the path in nighttime is higher than that for daytime.For example, a quasi-sinusoidal pattern is seen for NLK with nighttime amplitude of about 20 dB while daytime amplitude decreases by 10 dB.Other two stations have regular patterns with two peaks and two dips in the amplitude variation per day.The peaks with more fluctuations indicate the nighttime period and other peaks with smooth temporal dependence for daytime.The sharp dips are due to the interference of the different propagation modes and are so-called terminator times.Moreover, sharp drops in amplitude for long time period with amplitude reaching to −10 dB are due to the temporal shutdown of the VLF transmitters such as seen on 1, 2, 7 and 14 March in NPM transmitter and on 7, 12, and 15 March in NWC transmitter (i.e., not by the propagation effect in the EIWG (Earth-Ionosphere WaveGuide)).
VLF amplitude anomaly is clearly identifiable for the high-latitude path associated with the geomagnetic storms.The regular diurnal patterns break with small nighttime amplitude just after the storm onset.The small nighttime amplitude continues during the recovery phase.Other two paths (NPM: mid-latitude and NWC: low-mid-latitude) have little effect on VLF amplitude from the storms, however small decrease in nighttime amplitude can be seen in the mid-latitude path with a time delay (∼ 1 day in the nighttime amplitude (8 March) after the storm onset).
Figure 3a shows the first example of anomalies in VLF propagation due to ionospheric perturbations during a single (isolated in time) geomagnetic storm.The top panel indicates the daily values of Trend, while second and third panels show the Dispersion and NF respectively for the high-latitude path (NLK-KCH).These three parameters are calculated by using the signal amplitude with a time interval between 11:00 and 13:00 UT normalized with the corresponding standard deviation (σ ) before the current day (−15 to −1 day of the current day).Fourth and fifth panels indicate the geomagnetic conditions and are Dst and AE indices respectively.
As seen from Fig. 3a, significant VLF propagation anomaly is recognized.Trend and NF decrease significantly exceeding its 2σ value on 18 March.Dispersion is also enhanced remarkably exceeding their 2σ values on 17 March.On the other hand, Dst index drops suddenly on 17 March 2013 and reaches about −150 nT.AE index increases sharply when Dst drops which indicates the strong geomagnetic disturbance started from 17 March 2013 and continued for a few days.We consider the VLF anomalies observed around 17 March are due to the geomagnetic storm.
The 1-day time delay between the onset of the storm indicated by Dst and AE indices and Trend and NF is mainly due to the time resolution of the statistical parameter (1 day) with a limited daily time period of the data analysis (nighttime).Although the storm effect was slightly recognized by irregular diurnal pattern started during the time period of the analysis on the storm onset day (17 March) (not shown), this effect was not apparent in the trend and NF in Fig. 3a.In the case of this single event started on 17 March, Dst value started decreasing at around 08:00 UT, became −100 nT at around 16:00 UT and reached minimum value of −135 nT at around 20:00 UT.During the time period of data analysis for NLK (complete darkness in the transmitter-receiver path from 11:00 to 13:00 UT), the storm went on and the regular diurnal pattern was partially disturbed.The two statistical parameters such as Trend and NF using the mean value of the signal intensity did not significantly change but the dispersion increased significantly because variability is more sensitive to the change in the diurnal pattern.Therefore the storm effect was observed in Trend and NF during the next nighttime (18 March) during the recovery phase of the storm when the small signal amplitude was persistently observed.
Figure 3b shows the second example of anomalies in VLF propagation due to ionospheric perturbations during multiple geomagnetic storms (several storms closely situated in time) for high-latitude path (NLK-NSB).The time interval used for the data analysis is the same as Fig. 3a (from 11:00 to 13:00 UT) because of high-latitude propagation and same season.Trend and NF indicate the clear anomalies with two peaks exceeding 2σ corresponding to the days of each storm started on 7 and 9 March in 2012 according to sharp decreases of Dst index.Although Dispersion has only one enhancement exceeding 2σ on 9 March when the second storm starts, it was found that VLF anomalies are clearly recognizable in association with successive geomagnetic storms.Note that the second storm with a large change in Dst (< −100 nT) produced less change in trend (decrease in amplitude) rather than the first moderate storm (−100 nT < Dst < −50 nT).This is because of decreasing moving average amplitude due to the large drop of amplitude in the first storm.The real nighttime amplitude (complete darkness in the propagation path) on the day of the second storm (between 11:00 and 13:00 UT) from the diurnal pattern is much smaller than that during the first storm (Fig. 2).

Statistical results
In this subsection, statistical results of the VLF subionospheric anomalies during magnetic storms and a geomagnetically quiet time period are presented to elucidate the space weather effect in the lower ionospheric conditions.
The following criteria were used to identify the time period of anomalies during the geomagnetically active and quiet time periods based on Dst and AE indices.The geomagnetically active time period was defined from the days when Dst index gets smaller than −50 nT and its recovery around 0 nT.When multiple storms occurred close in time (within a few days), we consider two storms as one event and start and end days of storm period are the start day of the first storm and recovery day of the last storm respectively.On the other hand, geomagnetically quiet time period was chosen when AE index has a small value (smaller than about 130 nT in average) and remains almost constant.As a result, we examined eight and seven geomagnetic storm time periods for the years 2012 and 2013 respectively.The detailed information of 15 geomagnetic storms used in the data analysis was given in Table 2.In Table 2, the minimum value of Dst and the maximum value of AE during the storm time period were indicated along with the start date.On the other hand, seven and eight cases were chosen arbitrarily for the quiet time period in 2012 and 2013 respectively.
Based on the storm and geomagnetically quiet time periods selected in the previous paragraph, we calculated the occurrence rate of VLF anomalies for geomagnetically active and quite time periods.We examined two different conditions (A and B) for the statistical study to identify the VLF anomaly associated with a geomagnetic storm.First condition for the VLF anomaly is that all three VLF parameters (i.e., Trend, Dispersion, and NF) exhibit the anomaly at least once during the storm time period (i.e., timing from the storm onset to the end of recovery phase) and we define the statistical results based on this condition as case A. The second condition is similar to the case A but at least two of three VLF parameters indicate anomaly at least once in the storm time period, which is the case B.Here anomalies for different statistical parameters do not have to occur simultaneously (same day) but can occur on different days during the storms to satisfy the conditions.For both cases, the occurrence rate of anomalies for each transmitter station was calculated as the number of VLF path included anomalies divided by all available combination of the VLF paths.For example, if we have the amplitude data for a transmitter at five receiving stations for eight storms during the time period of analysis, we define the available number of paths as 40 by the number of receiving station times the number of storms.Then the occurrence rate of anomaly is 50 % if 20 from all available paths (40) indicate the VLF anomaly.
Figures 4 and 5 summarize the statistical results of the observed subionospheric anomalies during magnetic storms and rather quiet geomagnetic activity (arbitrary chosen) for two different cases A and B respectively.Figures 4a-c show the histograms of occurrence rate of VLF subionospheric anomalies for case A in relation with the geomagnetic storms for NLK, NPM and NWC transmitters respectively based on the data in 2012, while Fig. 4d-f shows the results for 2013 in the same format as Fig. 4a-c.
As seen from the Fig. 4, NLK (high-latitude paths) have the largest rate of anomalies (43-50 %) during the storm times for both years among three transmitter stations.Both NWC (north-south path) and NPM (mid-latitude) are found to have smaller occurrence rate than NLK as 19-31 % during storm time.On the other hand occurrence rate for no storm time period is ranging from 9 to 25 %, which is generally smaller than that for storm time period except NPM.The NPM has the occurrence rate of anomaly being comparable for storm and non-storm time periods for both years, whilst NLK has the largest difference in occurrence rate between storm and non-storm time period in 2012 more than 5 times.
Figures 5a-f show the histograms of occurrence rate of VLF anomalies for case B with the same format as Fig. 4.Although the occurrence rate remarkably increases for all cases, and general tendency of the statistical results does not vary significantly from the case A. For example, occurrence rate for the high-latitude path during storms increases to 63-71 % (case B) from 42-50 % (case A).However the rate of increasing the occurrence rate between the storm and nonstorm time periods is different.The occurrence rate for the non-storm time period (increasing 2-3 times than A) is superior to that for storm time period (increasing 1.5-2 times than A).

Discussion
In this paper we clearly show statistically the relationship between the space weather (i.e., geomagnetic storm) and VLF subionospheric anomaly.In particular, high-latitude transmitter-receiver path has a much larger probability to exhibit the anomaly rather than those of low latitude during the storm period indicated by a large negative drop signature in Dst index and enhancement of AE index.Since propagation distances between the transmitters and receivers studied here are not significantly different ranging from 6400 to 7600 km, the VLF anomaly can be mainly controlled by the length and intensity of ionospheric disturbance along the propagation path.Therefore NLK paths experiencing the largest distance in the high latitude part among three VLF transmitters lead to the largest observed occurrence rate of the anomaly due to the geomagnetic activity.
The occurrence of VLF anomalies for the quiet time period may indicate the anomalies other than space origin such as from atmosphere and lithosphere.Notably the occurrence rate for the non-storm time period is comparable to that of storm time period for the mid-latitude (NPM) and northsouth (NWC) paths, which suggest the geomagnetic activity gives little effects to the occurrence of the VLF anomaly.Moreover the mid-latitude path has the largest occurrence rate than others in quiet time period which indicates that anomalies other than space origins are significant at this latitude.
Severe meteorological process such as typhoons, cyclones, weather front, tornadoes and large thunderstorm systems may be the major sources of ionospheric anomalies in the atmosphere (troposphere) detected in our statistical parameters in the low-to mid-latitudes.These meteorological events are mostly distributed in the low-to mid-latitude and produce Atmospheric Gravity Waves (AGW) in the troposphere (Hung et al., 1979;Kelley, 1997;Laštovička, 2006;Xiao et al., 2007).The generated AGWs propagate upward and can disturb the ionosphere with a duration of hours (thunderstorm cells and tornadic storms) to even days (typhoon and cyclone).These long-lasting perturbation sources are suitable to be detected as ionospheric anomalies in the statistical parameters averaged over the nighttime period (Trend, Dispersion and NF).
Pre-seismic ionospheric anomaly is one of the lithospheric sources and is effectively identified as VLF anomaly in our statistical parameters because it may last for a relatively long time (at least for a night).The ionospheric anomaly of this type has been reported for years at various latitudes (lowto high-latitudes) (see e.g., Hayakawa and Hobara, 2010).These anomalies that occur a few days to 10 days prior to the earthquake are identified by the statistical parameters used in this study (Trend, Disersion and NF).Although preseismic anomalies could be included during the geomagnetically quiet time period, they may not be significant in this paper because VLF transmitter-receiver paths we used were mostly over the ocean (i.e., good conductor).The seawater above the epicenter may mask the pre-seismic electromagnetic phenomena and makes them difficult to propagate toward the sea surface from the epicenter then to the overlaying ionosphere.Therefore the transmitter-receiver paths over the continent can be more effective to detect the seismoionospheric anomalies.
In general many more VLF anomalies than those related to major seismic activities are observed considering the occurrence frequency.This means that the ionosphere is highly variable and various different perturbation sources exist.In this paper, we characterize the most pronounced effect (i.e., geomagnetic storm) to the VLF anomalies and will potentially use this result to separate the anomalies due to the geomagnetic storm from those from other unknown sources including pre-seismic origin during geomagnetically quiet days.Although it is still uncertain that only VLF data can be used to assess the probability of future seismic activities, the characteristics of VLF signals can form a part of the set of parameters that will allow us to provide a warning of possible hazards.Moreover two major physical mechanisms for such a coupling have been proposed and investigated such as elec-tric field (Pulinets et al., 2000) and AGW (Korepanov et al., 2009) but none of them have been experimentally proved yet.Therefore we need to investigate simultaneous ground and ionospheric data set to confirm if one of those/both mechanisms work to elucidate the observed ionospheric perturbations.
Other types of lithospheric short-lived sources such as earthquake main shock generates seismogenic AGW/TID (Fedorenko et al., 2005), associated Tsunami (Rozhnoi et al., 2012) as well as powerful explosion cause the ionospheric anomaly.These anomalies may occur at high latitude, however these anomalies do not last for a long time to affect the averaged nighttime properties of the ionosphere as we used.
Electromagnetic waves emitted from the tropospheric lightning propagate in the magnetosphere as plasma waves so-called lightning whistler.The lightning whistler causes the lightning-induced electron precipitations (LEP) due to the wave-particle interactions around the magnetic equator (e.g., Rycroft, 1973;Voss et al., 1984).These energetic particle precipitations cause the ionospheric perturbations in the D/E region and are detected as amplitude/phase changes of the VLF transmitter signals so-called classical Trimpi events (Helliwell et al., 1973) but each classical Trimpi event does not last for a very long time (e.g ∼ 100 s recovery time) and may not cause significant changes in our statistical parameters.Similar types of short-lived Trimpi events so-called early VLF events are observed with very small time delay from the causative lightning (<100 ms) in comparison to that of classical Trimpi (∼ 1 s) (Inan et al., 1988).These early VLF events are also sometimes associated with Transient luminous events (TLEs) such as red sprites and elves (Inan et al., 1995;Dowden et al., 1996;Hobara et al., 2001).A new class of early VLF event associated with lighting with long recovery time up to 20 min has been reported (Cotts and Inan, 2007).Even this long recovery event may not have a recovery time long enough to affect three nighttime averaged statistical parameters.
On the other hand, part of the ionosphere perturbations detected for the low and mid-latitude paths is due to the charge exchange mechanisms as well as wave-particle interactions in relation with magnetic storm, which may have a time period long enough to affect our statistical parameters (∼ longer than nighttime).The charge exchange mechanism suggests that the nighttime and even daytime E-region density enhances at low-mid-latitude due to the precipitated energetic neutral particles generated by the charge exchange of ring current ions during the storm activity (Lyons and Richmond, 1978;Voss and Smith, 1980;Jain and Singh, 1992).Precipitations of energetic electrons due to the pitch-angle scattering and diffusion due to the cyclotron resonance at the lower edge of the inner radiation belt may disturb D/E region ionosphere during the storm time period (Jain and Singh, 1990).
The occurrence rate significantly varies due to the conditions between the cases A and B. When we ease the condition from the strict case A (all three statistical parameters indicate anomaly during the storm period) to case B (at least two statistical parameters among three indicate anomaly during the storm time period), the occurrence rate for both storm and quiet times increases.However the amount of increase in the occurrence rate for the storm period (1.5 times) from the case A to case B is much smaller than that for the quiet time period (occurrence rate increases 2 times for many cases from the case A to B).This means that geomagnetic storms intrinsically affect VLF anomaly with all three statistical parameters, which are related each other.Therefore the occurrence rate does not increase so much even the condition is eased to case B from A. On the other hand, during the quiet time condition anomaly does not tend to occur for all three parameters (they are less correlated to each other).Therefore the occurrence rate increases much more when the condition is eased to case B from A.
Indeed during the geomagnetic storm time period, both VLF amplitude depression and the fluctuation were frequently observed.The observed amplitude decrease corresponds to the electron density enhancement D-region caused by high-energy auroral electron precipitation (Cummer et al., 1996(Cummer et al., , 1997)), which leads the change in Trend and NF values.Moreover amplitude fluctuations may indicate the movement of the auroral electrojet (Peter et al., 2006).Therefore our VLF parameters such as anomaly in Trend, NF corresponds to the decrease in amplitude depression, while increase in dispersion coincides with fluctuation.
The parameter "nighttime fluctuations" was defined differently in different research group such as described in Ray et al. (2011).In Ray et al. (2011), the nighttime fluctuations were defined as the significant deviations of signal amplitude from its mean value, which are calculated by using daily mean and standard deviations rather than running mean (previous 15 days) used in this paper.And the detected anomalies were correlated with major seismic activities near equatorial VLF paths.It is worthwhile analyzing our VLF data by using the nighttime fluctuations in Ray et al. (2011) in the future work to see any improvement in detection of ionospheric anomalies due to geomagnetic storm in comparison to this study.
Occurrence rate of anomaly during the storm time period between different years has similar tendency for NLK (high latitude) and NPM (mid-latitude).However the occurrence rate for NWC (north-south paths including low-midlatitudes) decrease significantly from 31 % in 2012 to 17 % in 2013 (2014) for the case A (smaller than the case for midlatitude path (NPM)).For the case B, this tendency is even more significant (from 62 % in 2012 to 22 % in 2013( 2014)).
The reason why such a big difference in the occurrence rate appeared is unknown because storm characteristics such as a max value of AE index and minimum value of Dst index shown in Table 2 are similar between studied years.Also average Trend value before the storm onset (averaged over a few days) does not vary significantly for all paths including NWC and NPM transmitters between 2012and 2013(2014).Hence the signal to noise ratio does not change significantly between the two time periods.Unlike high-latitude paths strongly affected by the auroral activity (these paths have always the largest detection rate for both years), lowmid-latitude range has much more complicated mechanism to get VLF anomalies in response to the energetic particle precipitations associated with geomagnetic storms by, e.g., charge exchange mechanism and wave-particle interactions in the magnetosphere.We continue investigating this point for future work.

Conclusions
We demonstrated the statistical properties of VLF signal amplitude based on 16-21 different paths consisting of seven receiving stations and three transmitters, which form highlatitude, mid-latitude and north-south paths for the time period of about 2 years.As a result, all three statistical parameters such as Trend, Dispersion and NF tend to indicate significant anomalies during the geomagnetic storm time period in the high-latitude paths (43-50 %) due to the auroral energetic electron precipitation.The mid-latitude and north-south paths have a little effect from the geomagnetic activity because of smaller occurrence rate of anomalies with similar occurrence rate as geomagnetically quiet time period (17-31 %).These anomalies are caused mainly by atmospheric and/or lithospheric origins.The obtained results are basic and important information not only to identify the space weather effects toward the lower ionosphere but also to separate various external physical causes of ionospheric disturbances.

Figure 1 .
Figure 1.Maps showing the geographic locations of VLF transmitters and receivers used for the data analysis.Corresponding great circle paths (GCPs) for every transmitter-receiver pair are also given.(a) Paths around Japan (indicating VLF receiving stations) and (b) entire paths (including VLF transmitters).

Figure 2 .
Figure 2. Diurnal patterns of VLF amplitude for high-latitude (NSB-NLK) (top panel), mid-latitude (NSB-NPM) (second panel) and northsouth path (NSB-NWC) (third panel) from 1 to 16 March in 2012.Local nighttime periods are shown in red color in the diurnal patterns.The bottom panel shows Dst variations for the relevant time period.

Figure 3 .
Figure 3. (a) Daily variations of VLF subionospheric propagation (three panels from the top), Dst (fourth panel) and AE indices (fifth panel) around a single (isolated in time) geomagnetic storm.Three statistical parameters for VLF propagation characteristics are trend (TR), dispersion (DP), and nighttime fluctuation (NF) from the top.(b) Same type of variation as (a) but for multiple geomagnetic storms (several storms occurred closely in time).

Figure 4 .
Figure 4. Histograms summarizing the statistical results for the case A. Every histogram contains the occurrence rate of VLF anomaly based on different statistical parameters for storm time and geomagnetically quiet time periods.(a) NLK, (b) NPM, (c) NWC for 2012, (d) NLK, (e) NPM, (f) NWC for 2013.

Figure 5 .
Figure 5. Histograms summarizing the statistical results for the case B with the same format as the case A in Fig. 4.

Table 1 .
Geophysical coordinates and transmitting frequencies of VLF/LF stations used in this paper, (a) transmitters and (b) receivers.

Table 2 .
Summary table indicating the geomagnetic storm events used in this paper.