Next Article in Journal
Leggett-Garg Inequalities for Quantum Fluctuating Work
Next Article in Special Issue
The Complexity Measures Associated with the Fluctuations of the Entropy in Natural Time before the Deadly México M8.2 Earthquake on 7 September 2017
Previous Article in Journal
Modulation Signal Recognition Based on Information Entropy and Ensemble Learning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Criticality Analysis of the Lower Ionosphere Perturbations Prior to the 2016 Kumamoto (Japan) Earthquakes as Based on VLF Electromagnetic Wave Propagation Data Observed at Multiple Stations

by
Stelios M. Potirakis
1,2,*,
Tomokazu Asano
2 and
Masashi Hayakawa
2
1
Department of Electrical and Electronics Engineering, University of West Attica, Campus 2, 250 Thivon and P. Ralli, Aigaleo, Athens, GR-12244, Greece
2
Hayakawa Institute of Seismo Electromagnetics Co. Ltd., University of Electro-Communications (UEC), Alliance Center #521, 1-1-1 Kojimacho, Chofu, Tokyo 182-0026, Japan
*
Author to whom correspondence should be addressed.
Entropy 2018, 20(3), 199; https://doi.org/10.3390/e20030199
Submission received: 17 February 2018 / Revised: 12 March 2018 / Accepted: 14 March 2018 / Published: 16 March 2018

Abstract

:
The perturbations of the ionosphere which are observed prior to significant earthquakes (EQs) have long been investigated and could be considered promising for short-term EQ prediction. One way to monitor ionospheric perturbations is by studying VLF/LF electromagnetic wave propagation through the lower ionosphere between specific transmitters and receivers. For this purpose, a network of eight receivers has been deployed throughout Japan which receive subionospheric signals from different transmitters located both in the same and other countries. In this study we analyze, in terms of the recently proposed natural time analysis, the data recorded by the above-mentioned network prior to the catastrophic 2016 Kumamoto fault-type EQs, which were as huge as the former 1995 Kobe EQ. These EQs occurred within a two-day period (14 April: M W = 6.2 and M W = 6.0 , 15 April: M W = 7.0 ) at shallow depths (~10 km), while their epicenters were adjacent. Our results show that lower ionospheric perturbations present critical dynamics from two weeks up to two days before the main shock occurrence. The results are compared to those by the conventional nighttime fluctuation method obtained for the same dataset and exhibit consistency. Finally, the temporal evolutions of criticality in ionospheric parameters and those in the lithosphere as seen from the ULF electromagnetic emissions are discussed in the context of the lithosphere-atmosphere-ionosphere coupling.

1. Introduction

In a large number of relevant articles published during the last few decades a variety of electromagnetic (EM) phenomena have been reported to appear prior to an earthquake (EQ), e.g., [1,2,3,4,5,6,7,8]. Specifically, the ionosphere has statistically been confirmed to be correlated with EQs in such a way that it presents remarkable sensitivity to EQ preparation processes happening in the lithosphere, e.g., [9,10]. A variety of pre-EQ EM phenomena related to different kinds of ionospheric anomalies have long been investigated at different frequency bands using multiple methods and could be considered promising for short-term EQ prediction [4,7,10,11,12,13,14,15,16,17,18,19,20].
The catastrophic 2016 Kumamoto EQs (14–15 April 2016) caused 50 deaths, over 1800 injuries and serious damage to local infrastructures [21]. These were fault-type EQs just like the 1995 Kobe EQ [22] and the magnitude of the main shock was as big as that of the Kobe EQ. In this sense, it is worthwhile to investigate the EM phenomena for these Kumamoto events. A variety of different ionosphere-related precursory signatures possibly related to the Kumamoto EQs have already been reported, such as ionospheric anomalies as detected by total electron content (TEC) calculated by GNSS receivers [23], or GIM-TEC [24], or GPS-TEC [25], transient variations on atmospheric and ionospheric parameters [26,27], subionospheric very low frequency (VLF) propagation anomalies [28,29], atmospheric ULF/ELF radiation and ULF depression [30], criticality in the ground-observed ULF magnetic fields [31], and intermittency-induced criticality in the VLF subionospheric propagation data [32].
We focus here on the lower ionosphere perturbations as monitored by the VLF propagation anomalies. A network of eight VLF/LF receivers has been operating during the last few years throughout Japan which receive subionospheric signals from different transmitters located both in the same and other countries. In this study we perform a criticality analysis of the data acquired by the specific network prior to the 2016 Kumamoto EQs. The analysis is performed by means of the recently proposed natural time (NT) method [33,34,35]. Beyond the seismic electric signals (SES) of Varotsos’ group, the NT method has already been successfully applied on other EM variations possibly related to EQs, such as MHz-kHz EM emissions [8,36,37] and ground-observed ULF magnetic fields (e.g., [31,38,39,40]). Yet, this is the first attempt of application of NT analysis to VLF subionospheric propagation data.
The way of application of the NT analysis is initially investigated, as well as its sensitivity to other possible causes of lower ionosphere anomalies. Then, the data from all eight stations of the aforementioned network of VLF/LF receivers concerning their receptions from a specific transmitter located in south Japan (with call sign “JJI”) are analyzed. The NT analysis results reveal that the lower ionosphere presented characteristics of critical dynamics from two weeks up to two days before the main shock. Moreover, as compared to the results obtained by analysis of the same data with the conventional nighttime fluctuation method and wave-hop theoretical computations [28,29] a remarkable consistency is observed.
The remaining of the article is organized as follows: Section 2 provides EQ information and describes the network of VLF/LF stations used in this study, as well as the data to be analyzed; a brief description of the key concepts and basic formulas of the NT analysis method is provided in Section 3; the NT analysis of the VLF subionospheric propagation data is presented in Section 4, along with a discussion on the sensitivity of the method to other causes of ionospheric anomalies. In Section 5, the results obtained in Section 4 are compared with the findings of the conventional nighttime fluctuation method, while the findings of NT analysis applied to the ULF electromagnetic emissions of the same period are discussed within the context of lithosphere-atmosphere-ionosphere (LAI) coupling; the global geomagnetic activity during the analyzed time period is also discussed. Finally, the main findings are summarized in Section 6.

2. EQs Information, VLF/LF Stations Network, Subionospheric Propagation Data

In this article we use VLF subionospheric propagation data between the transmitter “JJI”, located at Miyazaki in south Japan (geographic coordinates: 32.045° N, 130.811° E), shown in Figure 1 as a rectangle, and eight receivers dispersed all around Japan, shown in Figure 1 as triangles; the abbreviations of all receivers and all monitored transmitters along with their transmission signal frequencies are summarized in Table 1. Figure 1 also shows the 5th Fresnel zones only for three indicative links, defining an area over which the propagation path is considered sensitive to EQ preparation processes, e.g., [41]. In the same figure, all EQs with M W > 5.5 which happened in a wide area around Japan are also shown for the time period from 1 January 2016 to 30 April 2016. Even though an extremely low probability of ~1% was expected by the medium-term EQ prediction in the region of Kumamoto, a series of large magnitude EQs occurred, with a main event as big as the 1995 Kobe EQ. The three fault-type recent EQs of our interest happened in south-west Japan at a very close epicentral distance under the city of Kumamoto in the following order [42]: (1) M W = 6.2 , 14 April 2016, 12:26:41.1 UT (32.788° N, 130.704° E), depth ~ 9 km; (2) M W = 6.0 , 14 April 2016, 15:03:50.6 UT (32.697° N, 130.720° E), depth ~ 8 km; (3) M W = 7.0 , 15 April 2016, 16:25:15.7 UT (32.791° N, 130.754° E), depth ~ 10 km.
In this study we use the same data which are employed by the conventional nighttime fluctuation method [12]. Specifically, from the original raw measurements (sampling frequency f s = 1   Hz ) we first calculate the residue d A ( t ) between the received signal amplitude A ( t ) and an average signal amplitude A ( t ) calculated by means of a running average over ±15 days as d A ( t ) = A ( t ) A ( t ) . Then, since the (local) daytime data have been observed to exhibit too small disturbances to be analyzed, we use only (local) nighttime data over four time periods around the year (10:00–20:00 UT for 22/11–21/02, 11:00–19:00 UT for 22/02–21/05, 11:30–17:30 for 22/05–21/09, 10:30–19:00 for 22/09–21/11) to calculate daily values (1 value/day) for three quantities T R (“trend”), D P (“dispersion”), and N F (“nighttime fluctuation”):
T R = N s N e d A ( t ) N e N s ,
which is actually the mean value of d A ( t ) , where N s and N e are the time points of the start and end of the above defined nighttime periods,
D P = 1 N e N s N s N e ( d A ( t ) T R ) 2
which is actually the standard deviation of d A ( t ) , and
N F = N s N e ( d A ( t ) ) 2 .
Note that in the conventional nighttime fluctuation method, the normalized values of the above quantities, denoted respectively as D P * , T R * , and N F * , are usually studied. These are calculated as X * = ( X X ± 15 d a y s ) / σ ± 15 d a y s where X ± 15 d a y s and σ ± 15 d a y s denote the mean value and standard deviation ±15 days around the day of interest, respectively. However, in this paper we use the non-normalized values of D P , T R , and N F as defined in Equations (1)–(3). As an example of the observed VLF subionospheric propagation variations, Figure 2 shows TR* variation during the time period 1 January 2016–15 April 2016 for the receptions of the eight receivers from the transmitter JJI. We should clarify that the conventional nighttime fluctuation method is not employed in this paper; however, the already published conventional nighttime fluctuation method results related to the 2016 Kumamoto EQs [28,29] are discussed in the discussion section (see Section 5), in comparison with the herein presented results (see Section 4) which are obtained by the NT analysis method.

3. Natural Time Analysis Method

The natural time (NT) analysis method was originally proposed for the analysis of ultra-low frequency (≤1 Hz) SES signals [33,34,43], and has been shown to be optimal for enhancing the signals in the time-frequency space [44]. The transformation of a time series of “events” from the conventional time domain to the NT domain is performed by ignoring the time-stamp of each event and retaining only their normalized order (index) of occurrence. Explicitly, in a time series of N successive events, the NT, χ k , of the k th event is the index of occurrence of this event normalized by dividing by the total number of the considered events, χ k = k / N . On the other hand, the “energy”, Q k of each k th event is preserved. We note that the quantity Q k represents different physical quantities for various time series: for EQ time series it has been assigned to a seismic energy released (e.g., seismic moment) [34], and for SES signals that are of dichotomous nature it corresponds to SES pulse duration [34], while for geophysical scale MHz EM emission signals that are of non-dichotomous nature, it has been attributed to the energy of fracto-electromagnetic emission events as defined in [36]. The transformed time series ( χ k , Q k ) is then studied through the normalized power spectrum Π ( ϖ ) = | k = 1 N p k exp ( j ϖ χ k ) | 2 , where ϖ is the natural angular frequency, ϖ = 2 π φ , with φ standing for the frequency in NT, termed “natural frequency”, and p k = Q k / n = 1 N Q n corresponds to the k th event’s normalized energy. Note that the term “natural frequency” should not be confused with the rate at which a system oscillates when it is not driven by an external force; it defines an analysis domain dual to the NT domain, in the framework of Fourier–Stieltjes transform [35].
The study of Π ( ϖ ) at ϖ close to zero reveals the dynamic evolution of the time series under analysis. This is because all the moments of the distribution of p k can be estimated from Π ( ϖ ) at ϖ 0 [45]. Aiming to that, by the Taylor expansion Π ( ϖ ) = 1 κ 1 ϖ 2 + κ 2 ϖ 4 + the quantity κ 1 is defined, where κ 1 = k = 1 N p k χ k 2 ( k = 1 N p k χ k ) 2 , i.e., the variance of χ k weighted for p k characterizing the dispersion of the most significant events within the “rescaled” interval (0, 1]. Moreover, the entropy in NT, S n t , is defined [46] as S n t = k = 1 N p k χ k ln χ k ( k = 1 N p k χ k ) ln ( k = 1 N p k χ k ) and corresponds [35,46] to the value at q = 1 of the derivative of the fluctuation function f l ( q ) = χ q χ q with respect to q (while κ 1 corresponds to f l ( q ) for q = 2 ). It is a dynamic entropy depending on the sequential order of events [46]. The entropy, S n t , obtained upon considering [46] the time reversal T , i.e., T p m = p N m + 1 , is also taken into account.
A system is considered to approach criticality when the parameter κ 1 converges to the value κ 1 = 0.070 and at the same time both the entropy in NT and the entropy under time reversal satisfy the condition S n t , S n t < S u = ( ln 2 / 2 ) 1 / 4 [47], where S u stands for the entropy of a “uniform” distribution in NT [46]. It has to be mentioned that the criterion of the κ 1 = 0.070 value has originally been derived for SES activity and later on the basis of the Ising model. Its validity has been confirmed on real SES time series, while it has also been verified to be valid for several self-organized criticality (SOC) models and real time series of a variety of applications. In all these dynamical systems, it has been found that the value κ 1 = 0.070 can be considered as quantifying the extent of the organization of the system at the onset of the critical stage [35].
In the special case of NT analysis of foreshock seismicity [34,43,46,48], the seismicity is considered to be in a true critical state, a “true coincidence” is achieved, when three additional conditions are satisfied: (i) The “average” distance D between the curves of normalized power spectra Π ( ϖ ) of the evolving seismicity and the theoretical estimation of Π ( ϖ ) , Π c r i t i c a l ( ϖ ) = ( 18 / 5 ϖ 2 ) ( 6 cos ϖ / 5 ϖ 2 ) ( 12 sin ϖ / 5 ϖ 3 ) , Π c r i t i c a l ( ϖ ) 1 κ 1 ϖ 2 , for κ 1 = 0.070 should be smaller than 10 2 , i.e., D = | Π ( ϖ ) Π c r i t i c a l ( ϖ ) | < 10 2 (this is a practical criterion for signaling the achievement of spectral coincidence) [35]; (ii) the parameter κ 1 should approach the value κ 1 = 0.070 “by descending from above”, i.e., before the main event the parameter κ 1 should gradually decrease until it reaches the critical value 0.070 (this rule was found empirically) [35,43]; (iii) Since the underlying process is expected to be self-similar, the time of the true coincidence should not vary upon changing (within reasonable limits) either the magnitude threshold, M t h r e s , or the area, used in the calculation.
It should be finally clarified that in the case of seismicity analysis, the temporal evolution of the parameters κ 1 , S n t , S n t , and D is studied as new events that exceed the magnitude threshold M t h r e s are progressively included in the analysis. Specifically, as soon as one more event is included, first the time series ( χ k , Q k ) is rescaled in the NT domain, since each time the k th event corresponds to a NT χ k = k / N , where N is the progressively increasing (by each new event inclusion) total number of the considered successive events; then all the parameters involved in the NT analysis are calculated for this new time series; this process continues until the time of occurrence of the main event.
Note that in the case of NT analysis of foreshock seismicity, the introduction of magnitude threshold, M t h r e s , excludes some of the weaker EQ events (with magnitude below this threshold) from the NT analysis. On one hand, this is necessary in order to exclude events for which the recorded magnitude is not considered reliable; depending on the installed seismographic network characteristics, a specific magnitude threshold is usually defined to assure data completeness. On the other hand, the use of various magnitude thresholds, M t h r e s , offers a means of more accurate determination of the time when criticality is reached. In some cases, it happens that more than one time-points may satisfy the rest of NT critical state conditions, however the time of the true coincidence is finally selected by the last condition that “true coincidence should not vary upon changing (within reasonable limits) either the magnitude threshold, M t h r e s , or the area, used in the calculation”.

4. Analysis of Lower Ionosphere Nighttime Fluctuations Prior to the 2016 Kumamoto EQs

We apply here the NT method to the nighttime VLF propagation characteristic quantities T R , D P , and N F defined in Section 2, in a similar way to that for the ULF characteristics corresponding to magnetic field variation recorded prior to significant EQs [31,38,39,40]. Specifically: (i) We consider each daily value which is above a certain threshold as an event. In our nighttime VLF propagation characteristic quantities cases ( T R , D P , and N F ), the “energy” of k th event, that is the value of the quantity Q k of NT analysis (see Section 3), is considered to be equal to the corresponding non-normalized value of each one of the above quantities (as defined by Equations (1)–(3)), provided that this is above a certain threshold such as T R t h r e s , D P t h r e s , and N F t h r e s , respectively; (ii) Then, the NT analysis is performed as in the case of pre-EQ seismic activity on the revealed “events”. Starting from a specific day, all the parameters ( κ 1 , S n t , S n t , D defined in Section 3) are calculated for the time series of events rescaled in the NT domain each time a new event is added, checking for the corresponding criticality criteria as presented in Section 3 for the case of seismicity.
Before applying the NT analysis, the starting time of the analysis has to be determined. Since the ionosphere is known to be sensitive not only to pre-EQ processes, but also to a variety of different kinds of phenomena such as solar flares, magnetic storms, typhoons, tsunamis, and volcano eruptions, e.g., [41,49,50], we have first to check for the sensitivity of the NT analysis to such kind of phenomena. If the NT analysis is sensitive, i.e., if it indicates critical features before such phenomena, these should be either excluded from the analysis period, by setting the start time after their occurrence, or they have to be taken into account during the evaluation of the analysis results. This is because, if the NT analysis intended to examine the possibly EQ-related behavior of the ionosphere starts before a non-EQ-preparation-related phenomenon which undergoes a critical state, this might cause “masking” of the possible critical behavior of the ionosphere due to any EQ preparation processes.
Knowing that the conventional nighttime fluctuation method indicated a clear pre-EQ anomaly before the 2016 Kumamoto EQs in the JJI-IMZ path [28] (see Table 1 for station abbreviations), as well as the fact that a volcanic eruption occurred in the wider area (Sakurajima Volcano, geographic coordinates: 31°35′ N, 130°39′ E) on 5 February 2016 19:13 JST, we initially tried to apply NT analysis for the time period 1 January 2016–15 April 2016. The results obtained for the specific time period show that the NT analysis of the studied nighttime VLF propagation characteristics is affected by the volcanic eruption focusing on this phenomenon and masking any possible critical behavior due to phenomena which happened after that, including possible EQ-related ones. As an example, part of the NT analysis results for NF of the JJI-IMZ path during the time period 1 January 2016–15 April 2016 are shown in Figure 3, indicating that criticality is reached a few days before the abovementioned volcanic eruption. Specifically, according to the NT analysis criticality criteria for the case of seismicity (see Section 3) we observe in Figure 3 that criticality criteria are satisfied on 24 January 2016, i.e., ~2 weeks before the eruption. Therefore, since a disturbance in nighttime VLF propagation characteristics due to any phenomenon can be monitored by observing the fluctuation of their normalized version DP*, TR*, and NF* (see Section 2) it was decided to follow the rule of setting the initial time point for NT analysis at least a few (e.g., ~5) days after the day for which any normalized nighttime VLF propagation characteristic has exceeded the limit of ±2σ. According to this rule, it was decided that in order to reveal any possible critical characteristics related to the Kumamoto EQs, the NT analysis should start on 20 March 2016, because on 14 March 2016 there was high variation of the abovementioned normalized characteristics for some of the employed stations [28] (see also Figure 2).
The results of the NT analysis of N F , D P , and T R for all eight receiving stations and the time period 20 March 2016–15 April 2016 are summarized in Table 2, while indicative results are shown in Figure 4, Figure 5, Figure 6 and Figure 7. In Figure 4, Figure 5 and Figure 6 clear criticality indications for N F , D P , and T R , respectively, are presented, while Figure 7 portrays examples of marginal criticality indications (please refer to NT analysis criticality criteria for the case of seismicity in Section 3). For example, we can see in Figure 4 that there is a specific time period, namely 11–12 April 2016, during which the NT analysis criticality criteria, i.e., κ 1 reaching the value 0.070 “from above” while at the same time S n t , S n t < ( ln 2 / 2 ) 1 / 4   ( 0.0966 ) and D < 10 2 are satisfied for the presented four different threshold cases. Therefore, the analysis of N F for the JJI-IMZ path during the time period 20 March 2016–15 April 2016 which is presented in Figure 4 shows clearly that the specific quantity reaches criticality on 11–12 April 2016. As an example of marginal criticality indications, we can refer to the analysis of D P for the JJI-ANA path shown in Figure 7d. In the specific case, we consider that the criterion that the parameter κ 1 should approach the value κ 1 = 0.070 “by descending from above” is marginally satisfied, since only a few values of κ 1 could be calculated for the specific threshold before it approached the value κ 1 = 0.070 on 28–30 March 2016. We should clarify that the existence of just one threshold (or very few thresholds) for which criticality conditions are satisfied, and NT analysis parameters ( κ 1 , S n t , S n t , D ) values which are very close to the limits set by the corresponding criteria, are also considered as “marginal criticality indications”. We observe from Table 2 as follows:
(1)
Criticality has been revealed for all eight stations, but for different propagation characteristics ( N F , D P , and T R ) we have criticality for different combinations of stations and dates.
(2)
It is especially worth noting that the receiving stations KMK, TYH, ANA and KTU, which are all situated on the east (Pacific) coast of Japan showed either marginal or clear indications of criticality from 28 March 2016 up to 1–2 April 2016. However, it may be possible that this behavior is related to the M5.9 EQ which happened in the Pacific Ocean coast of Japan on 1 April 2016, 02:39:08.06 UT (33.3835° N, 136.3857° E), depth = 14 km (please also see Figure 1), and not to the 2016 Kumamoto EQs.
(3)
After those dates, clear criticality indications are progressively appearing from 5 April 2016 (i.e., 9–10 days before the 2016 Kumamoto EQs) up to 13 April 2016 (i.e., 1–2 days before the 2016 Kumamoto EQs), while marginal indications of criticality are observed even on 15 April 2016 but only at KMK station.
(4)
The timeline of clear criticality indications develops as follows: (a) first criticality is found in D P of the JJI-KTU path (5–7 April 2016); (b) then, criticality appears in the JJI-STU path (starting from 8 April 2016 in D P , and T R ); (c) between 9 April 2016 and 12–13 April 2016, clear criticality evidence appear in five propagation paths of JJI-STU, JJI-KTU, JJI-IMZ, JJI-KMK, JJI-AKT, some of them presenting critical characteristics in more than one of the analyzed quantities ( T R , D P , and N F ).

5. Discussion

First we try to compare the NT analysis results presented in the previous section with the results by the conventional analysis presented in [28]. As indicated by points (2) and (3) in the previous section, we pay attention to the time period after 5 April because propagation anomalies during this period are more likely to be related with the Kumamoto EQs. In [28] it has been reported that propagation anomalies were observed during the time period 19 March–10 April, which is consistent with the NT analysis results presented in this paper.
A further extensive study has recently been presented in [29]. The temporal evolutions of amplitude changes at all eight receiving stations, have enabled Asano and Hayakawa [29] to compare those with the theoretical estimations with the use of wave-hop method. In the wave-hop method, one can change, independently, the reflection height (lower ionosphere height) (either increasing or decreasing) of 1-hop sky wave and that of 2-hop sky wave close to the transmitter, to estimate the amplitude at all stations. The comparisons between the observed and theoretical amplitudes at all stations allowed them to deduce the spatio-temporal evolution of the ionospheric perturbation associated with the Kumamoto EQs. As a result, the perturbation begins to appear on 3–4 April 2016 (about two weeks before the EQs) and it continues to develop spatially, i.e., horizontal spatial extent is expanding and vertical scale is decreasing (decrease in the lower ionosphere). The maximum of spatial development happens 10–12 April 2016, and is followed by a rapid decay. These spatio-temporal evolutions are found to be in extremely good consistency with the NT analysis findings of criticality in this paper. The dates of maximum spatial development in the ionospheric perturbation are consistent with those in Table 2 indicating when criticality was reached for different propagation parameters of different subionospheric paths.
The generation mechanism of seismo-ionospheric perturbation is poorly understood, but a few hypotheses have already been proposed [50,51,52,53]. The major agent of ionospheric perturbations is likely to be located in the lithosphere or in the near surface. So, it is worthwhile to compare the NT analysis results of this paper with those referring to the lithosphere as seen from the lithospheric ULF electromagnetic emissions as recorded by the ground-based magnetic observatory of Kanoya (KNY) [31] (beyond the time series analysis, full information about the observatory and the involved instrumentation is also provided in [31]). Note that in both cases NT analysis was applied to time series of daily values (1 value/day), representing VLF subionospheric propagation characteristics and ULF magnetic field characteristics, respectively. Figure 8 illustrates a comparison of dates of criticality revealed by the NT analysis method in the ionosphere (top panel) and that in the lithosphere as seen from ULF electromagnetic radiation (bottom panel). This figure indicates that criticality in the lithosphere has been reached about 1 month to 2 weeks before the Kumamoto EQs, whereas criticality can be observed in the ionosphere as VLF propagation anomaly about 2 weeks to a few days before the EQs. So we can find a significant difference in dates of criticality (of the order of 1–2 weeks) in the lithosphere and in the ionosphere. This precedence concerning the appearance of critical dynamics in the ULF magnetic field variations as compared to the appearance of criticality in the lower ionosphere could imply that the mechanism producing the ULF magnetic field anomaly drives the mechanism producing the VLF subionospheric propagation anomaly. This might be of great importance in studying the mechanism of LAI coupling. However, a further investigation of this issue is considered out of the scope of this paper and will be pursued in the future.
As it is well known, magnetosphere’s condition influences ionosphere. Geomagnetic disturbances such as storms, sudden commencements etc. highly influence ionosphere. On the other hand, a new interpretation on the relation between K p index and EQs has been suggested in recent studies [54,55]. The planetary K p index is obtained from a number of magnetometer stations at mid-latitudes and reflects global geomagnetic activity. Moreover, the main indices employed for the monitoring of magnetic effects, K p , D s t , and A E , correlate rather well during periods of noticeable disturbances. During the analyzed period, 20 March 2016–15 April 2016, there were no noticeable magnetic field disturbances (max three-hourly K p reached the value 6- on 7 April 2016), although there were variations with relative increase of K p during specific periods as well as sudden commencements (http://www.gfz-potsdam.de/en/kp-index). It is not easy to attempt any correlation of the obtained NT analysis results with the variation of K p , however it would be very interesting a direct application of the NT analysis to the K p data and a comparison with the ionospheric NT analysis results in a future study.

6. Conclusions

The results of the first attempt of applying the NT criticality analysis to the subionospheric VLF data have been presented in the present article. As summarized in Section 4, the lower ionosphere as seen by VLF propagation has exhibited critical characteristics from two weeks up to two days before the main shock of the disastrous 2016 Kumamoto EQs (15 April 2016). We note that four of the receiving stations, which are all situated on the east (Pacific) coast of Japan, showed either marginal or clear indications of criticality from 28 March 2016 up to 1–2 April 2016. However, it may be possible that this behavior is related to the M5.9 EQ which happened in the Pacific Ocean coast of Japan on 1 April 2016, and not to the 2016 Kumamoto EQs.
Importantly, the NT analysis findings of criticality in this paper are found to be in extremely good consistency with the results of the conventional nighttime fluctuation method obtained for the same dataset and the theoretical calculations by means of the wave-hop method. Specifically, the time period for which VLF subionospheric propagation anomalies were identified by the conventional nighttime fluctuation method overlaps with the criticality dates revealed by the NT analysis method, while the spatio-temporal evolution of the ionospheric perturbation associated with the Kumamoto EQs obtained by the wave-hop method matches the progressive appearance of critical dynamics in the studied receivers.
We have also observed that the appearance of critical dynamics in the ground-observed ULF magnetic field recorded at a magnetic observatory close to the Kumamoto EQs epicenters precedes the appearance of criticality in the lower ionosphere. This could imply that the mechanism producing the ULF magnetic field anomaly drives the mechanism producing the VLF subionospheric propagation anomaly. However, more investigation is necessary. A multi-parameter analysis including the ULF/ELF radiation as a signature of atmospheric perturbation [56] is highly required in order to elucidate the mechanism of LAI coupling in the future. Moreover, it would be very interesting a direct application of the NT analysis to K p data and a comparison with the ionospheric NT analysis results in a future study.

Acknowledgments

Two of the authors (S.M.P. & M.H.) would like to express their gratitude to the Matsumae International Foundation for partly supporting this research by providing S.M.P. a 4-month research fellowship [Gno: 17G06] for visiting Hayakawa Institute of Seismo Electromagnetics Co. Ltd. Also thanks are due to Earthquake Analysis Laboratory for permitting us to use the VLF data as well as to the Northern California Earthquake Data Center and the Advanced National Seismic System (ANSS) for the access to the ANSS Composite Earthquake Catalog.

Author Contributions

T.A. and M.H. acquired the experimental data; S.M.P. and T.A. analyzed the data; S.M.P. and M.H. evaluated the analysis results; S.M.P. and M.H. wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

References

  1. Hayakawa, M.; Molchanov, O.A. (Eds.) Seismo Electromagnetics: Lithosphere-Atmosphere-Ionosphere Coupling; TERRAPUB: Tokyo, Japan, 2002; ISBN 4-88704-130-6. [Google Scholar]
  2. Pulinets, S.; Boyarchuk, K. Ionospheric Precursors of Earthquakes; Springer: Berlin, Germany, 2004; ISBN 978-3-540-26468-2. [Google Scholar]
  3. Varotsos, P.A. The Physics of Seismic Electric Signals; TERRAPUB: Tokyo, Japan, 2005; ISBN 4-88704-136-5. [Google Scholar]
  4. Molchanov, O.A.; Hayakawa, M. Seismo Electromagnetics and Related Phenomena: History and Latest Results; TERRAPUB: Tokyo, Japan, 2008; ISBN 978-4-88704-143-1. [Google Scholar]
  5. Hattori, K. ULF geomagnetic changes associated with major earthquakes. In Earthquake Prediction Studies: Seismo Elecromagnetics; Hayakawa, M., Ed.; TERRAPUB: Tokyo, Japan, 2013; pp. 129–152. ISBN 978-4-88704-163-9. [Google Scholar]
  6. Eftaxias, K.; Potirakis, S.M. Current challenges for pre-earthquake electromagnetic emissions: Shedding light from micro-scale plastic flow, granular packings, phase transitions and self-affinity notion of fracture process. Nonlinear Process. Geophys. 2013, 20, 771–792. [Google Scholar] [CrossRef]
  7. Hayakawa, M. Earthquake Prediction with Radio Techniques; John Willey and Sons: Singapore, 2015; ISBN 978-1-118-77016-0. [Google Scholar]
  8. Potirakis, S.M.; Contoyiannis, Y.; Melis, N.S.; Kopanas, J.; Antonopoulos, G.; Balasis, G.; Kontoes, C.; Nomicos, C.; Eftaxias, K. Recent seismic activity at Cephalonia (Greece): A study through candidate electromagnetic precursors in terms of non-linear dynamics. Nonlinear Process. Geophys. 2016, 23, 223–240. [Google Scholar] [CrossRef]
  9. Liu, J.Y. Earthquake precursors observed in the ionospheric F-region. In Electromagnetic Phenomena Associated with Earthquakes; Hayakawa, M., Ed.; Transworld Research Network: Trivandrum, India, 2009; pp. 187–204. ISBN 978-81-7895-297-0. [Google Scholar]
  10. Hayakawa, M.; Kasahara, Y.; Nakamura, T.; Muto, F.; Horie, T.; Maekawa, S.; Hobara, Y.; Rozhnoi, A.; Solovieva, M.; Molchanov, O.A. A statistical study on the correlation between lower ionospheric perturbations as seen by subionospheric VLF/LF propagation and earthquakes. J. Geophys. Res. Space Phys. 2010, 115. [Google Scholar] [CrossRef]
  11. Kon, S.; Nishihashi, M.; Hattori, K. Ionospheric anomalies possibly associated with M ≥ 6.0 earthquakes in the Japan area during 1998–2010: Case studies and statistical study. J. Asian Earth Sci. 2011, 41, 410–420. [Google Scholar] [CrossRef]
  12. Hayakawa, M. Probing the lower ionospheric perturbations associated with earthquakes by means of subionospheric VLF/LF propagation. Earthq. Sci. 2011, 24, 609–637. [Google Scholar] [CrossRef]
  13. Chmyrev, V.M.; Isaev, N.V.; Bilichenko, S.V.; Stanev, G.A. Observation by space-borne detectors of electric fields and hydromagnetic waves in the ionosphere over on earthquake center. Phys. Earth Planet. Int. 1989, 57, 110–114. [Google Scholar] [CrossRef]
  14. Fraser-Smith, A.C.; Bernardi, A.; McGill, P.R.; Ladd, M.E.; Helliwell, R.A.; Viilard, O.G. Low-frequency magnetic field measurements near the epicenter of the Ms 7.1 Loma Prieta Earthquake. Geophys. Rev. Lett. 1990, 17, 1465–1468. [Google Scholar] [CrossRef]
  15. Kopytenko, Y.A.; Matiashvili, T.G.; Voronov, P.M.; Kopytenko, E.A.; Molchanov, O.A. Detection of ultra-low frequency emissions connected with the Spitak earthquake and its aftershock activity, based on geomagnetic pulsations data at Dusheti and Vardzia observatories. Phys. Earth Planet. Int. 1993, 77, 85–95. [Google Scholar] [CrossRef]
  16. Aleksandrin, S.Y.; Galper, A.M.; Grishantzeva, L.A.; Koldashov, S.V.; Maslennikov, L.V.; Murashov, A.M.; Picozza, P.; Sgrigna, V.; Voronov, S.A. High-energy charged particle bursts in near-Earth space as earthquake precursors. Ann. Geophys. 2003, 21, 597–602. [Google Scholar] [CrossRef]
  17. Parrot, M. (Ed.) First results of the DEMETER micro-satellite. Planet. Space Sci. 2006, 54, 411–558. [Google Scholar]
  18. Liu, J.Y.; Chen, C.H.; Chen, Y.I.; Yang, W.H.; Oyama, K.I.; Kuo, K.W. A statistical study of ionospheric earthquake precursors monitored by using equatorial ionization anomaly of GPS TEC in Taiwan during 2001–2007. J. Asian Earth Sci. 2010, 39, 76–80. [Google Scholar] [CrossRef]
  19. Anagnostopoulos, G.C.; Vassiliadis, E.; Pulinets, S. Characteristics of flux-time profiles, temporal evolution, and spatial distribution of radiation-belt electron precipitation bursts in the upper ionosphere before great and giant earthquakes. Ann. Geophys. 2012, 55, 21–36. [Google Scholar] [CrossRef]
  20. Parrot, M. Satellite observations of ionospheric perturbations related to seismic activity. In Earthquake Prediction Studies: Seismo Electromagnetics; Hayakawa, M., Ed.; TERRAPUB: Tokyo, Japan, 2013; pp. 1–16. ISBN 978-4-88704-163-9. [Google Scholar]
  21. Sano, Y.; Takahata, N.; Kagoshima, T.; Shibata, T.; Onoue, T.; Zhao, D. Groundwater helium anomaly reflects strain change during the 2016 Kumamoto earthquake in Southwest Japan. Nat. Sci. Rep. 2016, 6, 37939. [Google Scholar] [CrossRef] [PubMed]
  22. Nagao, T.; Enomoto, Y.; Fujinawa, Y.; Hata, M.; Hayakawa, M.; Huang, Q.; Izutsu, J.; Kushida, Y.; Maeda, K.; Oike, K.; et al. Electromagnetic anomalies associated with 1995 Kobe earthquake. J. Geodyn. 2002, 33, 401–411. [Google Scholar] [CrossRef]
  23. Iwata, T.; Umeno, K. Preseismic ionospheric anomalies detected before the 2016 Kumamoto earthquake. J. Geophys. Res. 2017, 122, 3602–3616. [Google Scholar] [CrossRef]
  24. Kobari, T.; Han, P.; Hattori, K. b-value and GIM-TEC variations for Tokachi region and the 2016 Kumamoto earthquakes. In Proceedings of the International Workshop on Earthquake Preparation Process 2016—Observation, Validation, Modeling, Forecasting (IWEP2016), Chiba University, Chiba, Japan, 27–28 May 2016. [Google Scholar]
  25. Yagmur, M.; Hirooka, S.; Hattori, K. 3D Structure of ionospheric disturbances related to large earthquakes. In Proceedings of the International Workshop on Earthquake Preparation Process 2016—Observation, Validation, Modeling, Forecasting (IWEP2016), Chiba University, Chiba, Japan, 27–28 May 2016. [Google Scholar]
  26. Ouzounov, D.; Hattori, K.; Pulinets, S.; Petrov, L. Observation of transient signatures in atmosphere and ionosphere prior to the 2016 Kumamoto earthquake in Japan. Preliminary results. In Proceedings of the Japan Geoscience Union Meeting 2016, Chiba, Japan, 22–26 May 2016. [Google Scholar]
  27. Ouzounov, D.; Pulinets, S.; Hattori, K.; Liu, T.; Kafatos, M. An interdisciplinary approach of utilizing pre-earthquake signals for short-term prediction studies. In Proceedings of the 35th General Assembly of the European Seismological Commission, Trieste, Italy, 4–10 September 2016. ESC2016-401-1. [Google Scholar]
  28. Hayakawa, M.; Asano, T. Subionospheric VLF propagation anomaly prior to the Kumamoto earthquake in April, 2016. New Concepts Glob. Tecton. J. 2016, 4, 273–276. [Google Scholar]
  29. Asano, T.; Hayakawa, M. On the lower ionospheric perturbation for the 2016 Kumamoto earthquakes on the basis of VLF propagation data observed at multiple stations and wave-hop theoretical computations. In Proceedings of the International Workshop on Earthquake Preparation Process 2017—Observation, Validation, Modeling, Forecasting (IWEP2017), Chiba University, Chiba, Japan, 26–27 May 2017. [Google Scholar]
  30. Schekotov, A.; Izutsu, J.; Asamo, T.; Potirakis, S.M.; Hayakawa, M. Electromagnetic precursors to the 2016 Kumamoto earthquakes. Open J. Earthq. Res. 2017, 6, 168–179. [Google Scholar] [CrossRef]
  31. Potirakis, S.M.; Schekotov, A.; Asano, T.; Hayakawa, M. Natural time analysis on the ultra-low frequency magnetic field variations prior to the 2016 Kumamoto (Japan) earthquakes. J. Asian Earth Sci. 2018, 154, 419–427. [Google Scholar] [CrossRef]
  32. Potirakis, S.M.; Contoyiannis, Y.; Asano, T.; Hayakawa, M. Intermittency-induced criticality in the lower ionosphere prior to the 2016 Kumamoto earthquakes as embedded in the VLF propagation data observed at multiple stations. Tectonophysics 2018, 722, 422–431. [Google Scholar] [CrossRef]
  33. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S. Long-range correlations in the electric signals that precede rupture. Phys. Rev. E 2002, 66, 011902. [Google Scholar] [CrossRef] [PubMed]
  34. Varotsos, P.A.; Sarlis, N.V.; Tanaka, H.K.; Skordas, E.S. Similarity of fluctuations in correlated systems: The case of seismicity. Phys. Rev. E 2005, 72, 041103. [Google Scholar] [CrossRef] [PubMed]
  35. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S. Natural Time Analysis: The New View of Time; Springer: Berlin, Germany, 2011; ISBN 978-3-642-16449-1. [Google Scholar]
  36. Potirakis, S.M.; Karadimitrakis, A.; Eftaxias, K. Natural time analysis of critical phenomena: The case of pre-fracture electromagnetic emissions. Chaos Interdiscip. J. Nonlinear Sci. 2013, 23, 023117. [Google Scholar] [CrossRef] [PubMed]
  37. Potirakis, S.M.; Contoyiannis, Y.; Eftaxias, K.; Koulouras, G.; Nomicos, C. Recent field observations indicating an earth system in critical condition before the occurrence of a significant earthquake. IEEE Geosci. Remote Sens. Lett. 2015, 12, 631–635. [Google Scholar] [CrossRef]
  38. Hayakawa, M.; Schekotov, A.; Potirakis, S.M.; Eftaxias, K.; Li, Q.; Asano, T. An integrated study of ULF magnetic field variations in association with the 2008 Sichuan earthquake, on the basis of statistical and critical analyses. Open J. Earthq. Res. 2015, 4, 85–93. [Google Scholar] [CrossRef]
  39. Hayakawa, M.; Schekotov, A.; Potirakis, S.; Eftaxias, K. Criticality features in ULF magnetic fields prior to the 2011 Tohoku earthquake. Proc. Jpn. Acad. Ser. B 2015, 91, 25–30. [Google Scholar] [CrossRef] [PubMed]
  40. Potirakis, S.M.; Eftaxias, K.; Schekotov, A.; Yamaguchi, H.; Hayakawa, H. Criticality features in ultra-low frequency magnetic fields prior to the 2013 M6.3 Kobe earthquake. Ann. Geophys. 2016, 59, 0317. [Google Scholar] [CrossRef]
  41. Rozhnoi, A.; Solovieva, M.; Hayakawa, M. VLF/LF signals method for searching of electromagnetic earthquake precursors. In Earthquake Prediction Studies: Seismo Electromagnetics; Hayakawa, M., Ed.; TERRAPUB: Tokyo, Japan, 2013; pp. 31–48. ISBN 978-4-88704-163-9. [Google Scholar]
  42. Himematsu, Y.; Furuya, M. Fault source model for the 2016 Kumamoto earthquake sequence based on ALOS-2/PALSAR-2 pixel-offset data: Evidence for dynamic slip partitioning. Earth Planets Space 2016, 68, 169. [Google Scholar] [CrossRef]
  43. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S. Spatio-temporal complexity aspects on the interrelation between seismic electric signals and seismicity. Practica Athens Acad. 2001, 76, 294–321. [Google Scholar]
  44. Abe, S.; Sarlis, N.V.; Skordas, E.S.; Tanaka, H.K.; Varotsos, P.A. Origin of the usefulness of the natural-time representation of complex time series. Phys. Rev. Lett. 2005, 94, 170601. [Google Scholar] [CrossRef] [PubMed]
  45. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S.; Uyeda, S.; Kamogawa, M. Natural time analysis of critical phenomena. Proc. Natl. Acad. Sci. USA 2011, 108, 11361–11364. [Google Scholar] [CrossRef]
  46. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S.; Tanaka, H.K.; Lazaridou, M.S. Entropy of seismic electric signals: Analysis in the natural time under time reversal. Phys. Rev. E 2006, 73, 031114. [Google Scholar] [CrossRef] [PubMed]
  47. Sarlis, N.V.; Skordas, E.S.; Varotsos, P.A. Similarity of fluctuations in systems exhibiting self-organized criticality. Europhys. Lett. 2011, 96, 28006. [Google Scholar] [CrossRef]
  48. Sarlis, N.V.; Skordas, E.S.; Lazaridou, M.S.; Varotsos, P.A. Investigation of seismicity after the initiation of a Seismic Electric Signal activity until the main shock. Proc. Jpn. Acad. Ser. B 2008, 84, 331–343. [Google Scholar] [CrossRef]
  49. Rozhnoi, A.; Solovieva, M.; Levin, B.; Hayakawa, M.; Fedun, V. Meteorological effects in the lower ionosphere as based on VLF/LF signal observations. Nat. Hazards Earth Syst. Sci. 2014, 14, 2671–2679. [Google Scholar] [CrossRef]
  50. Hayakawa, M.; Asano, T.; Rozhnoi, A.; Solovieva, M. VLF/LF sounding of ionospheric perturbations and possible association with earthquakes. In Pre-Earthquake Processes: A Multi-Disciplinary Approach to Earthquake Prediction Studies; Ouzounov, D., Pulinets, S., Hattori, K., Taylor, P., Eds.; Willey: Indianapolis, IN, USA, 2018; ISBN 978-1-119-15694-9. [Google Scholar]
  51. Korepanov, V.; Hayakawa, M.; Yampolski, Y.; Lizunov, G. AGW as a seismo-ionospheric coupling responsible agent. Phys. Chem. Earth Parts A/B/C 2009, 34, 485–495. [Google Scholar] [CrossRef]
  52. Freund, F. Toward a unified solid state theory for pre-earthquake signals. Acta Geophys. 2010, 58, 719–766. [Google Scholar] [CrossRef]
  53. Pulinets, S.; Ouzounov, D. Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) model—An unified concept for earthquake precursors validation. J. Asian Earth Sci. 2011, 41, 371–382. [Google Scholar] [CrossRef]
  54. Anagnostopoulos, G.C.; Papandreou, A. Space conditions during a month of a sequence of six M > 6.8 earthquakes ending with the tsunami of 26 December 2004. Nat. Hazards Earth Syst. Sci. 2012, 12, 1551–1559. [Google Scholar] [CrossRef]
  55. Urata, N.; Duma, G.; Freund, F. Geomagnetic Kp Index and Earthquakes. Open J. Earthq. Res. 2018, 7, 39–52. [Google Scholar] [CrossRef]
  56. Schekotov, A.; Hayakawa, M. ULF/ELF Electromagnetic Phenomena for Short-term Earthquake Prediction; LAP LAMBERT Academic Publishing: Beau Bassin, Mauritius, 2017; ISBN 978-3-330-06286-3. [Google Scholar]
Figure 1. Map of the wider area of Japan, showing the network of VLF/LF receivers (triangles), the VLF transmitters (rectangles), the 5th Fresnel zones for the JJI-NSB, JJI-STU and JJI-KTU paths, as well as all EQs with M W > 5.5 which happened during the time period from 1 January 2016 to 30 April 2016. The circle size is proportional to EQ’s magnitude and its color refers to hypocenter’s depth. The date of occurrence appears only for the EQs with hypocenter at depth <100 km. The Kumamoto EQs correspond to the large-sized overlapping circles in south-west Japan (on Kyushu Island).
Figure 1. Map of the wider area of Japan, showing the network of VLF/LF receivers (triangles), the VLF transmitters (rectangles), the 5th Fresnel zones for the JJI-NSB, JJI-STU and JJI-KTU paths, as well as all EQs with M W > 5.5 which happened during the time period from 1 January 2016 to 30 April 2016. The circle size is proportional to EQ’s magnitude and its color refers to hypocenter’s depth. The date of occurrence appears only for the EQs with hypocenter at depth <100 km. The Kumamoto EQs correspond to the large-sized overlapping circles in south-west Japan (on Kyushu Island).
Entropy 20 00199 g001
Figure 2. VLF subionospheric propagation characteristic quantity TR* during the time period 1 January 2016–15 April 2016 for the receptions of the eight receivers from the transmitter JJI (see also Table 1 and Figure 1). Values exceeding −2σ threshold (shown as red horizontal lines) are considered as statistically significant anomalies (grey-color filled areas) according to the conventional nighttime fluctuation method [12]. Black vertical patches indicate time periods for which data are missing due to any reason. The horizontal time (x-) axis shows date (UT).
Figure 2. VLF subionospheric propagation characteristic quantity TR* during the time period 1 January 2016–15 April 2016 for the receptions of the eight receivers from the transmitter JJI (see also Table 1 and Figure 1). Values exceeding −2σ threshold (shown as red horizontal lines) are considered as statistically significant anomalies (grey-color filled areas) according to the conventional nighttime fluctuation method [12]. Black vertical patches indicate time periods for which data are missing due to any reason. The horizontal time (x-) axis shows date (UT).
Entropy 20 00199 g002
Figure 3. NT analysis of N F for the JJI-IMZ path for the time period 1 January 2016–15 April 2016. Variations of the NT analysis parameters ( κ 1 , S n t , S n t , and D ) for the different thresholds (a) 2, (b) 8; (c) 12; and (d) 20, respectively. The entropy limit of S u ( 0.0966 ) , the κ 1 value 0.070 and a region of ±0.005 around it are shown by the horizontal solid light green, solid grey and the grey dashed lines, respectively, while the D limit (10−2) is shown by the horizontal solid brown line. Note that the events employed depend on the considered threshold. Moreover, the time (x-) axis (date, UT) is not linear in terms of the conventional time of occurrence of the events, since the employed events appear equally spaced relative to x-axis as the NT representation demands, although they are not equally spaced in conventional time.
Figure 3. NT analysis of N F for the JJI-IMZ path for the time period 1 January 2016–15 April 2016. Variations of the NT analysis parameters ( κ 1 , S n t , S n t , and D ) for the different thresholds (a) 2, (b) 8; (c) 12; and (d) 20, respectively. The entropy limit of S u ( 0.0966 ) , the κ 1 value 0.070 and a region of ±0.005 around it are shown by the horizontal solid light green, solid grey and the grey dashed lines, respectively, while the D limit (10−2) is shown by the horizontal solid brown line. Note that the events employed depend on the considered threshold. Moreover, the time (x-) axis (date, UT) is not linear in terms of the conventional time of occurrence of the events, since the employed events appear equally spaced relative to x-axis as the NT representation demands, although they are not equally spaced in conventional time.
Entropy 20 00199 g003
Figure 4. Variations of the NT analysis parameters ( κ 1 , S n t , S n t and D ) corresponding to the analysis of NF for the JJI-IMZ path during the time period 20 March 2016–15 April 2016 for four threshold values (a) 4, (b) 7, (c) 11 and (d) 14, respectively. Figure format is similar to Figure 3.
Figure 4. Variations of the NT analysis parameters ( κ 1 , S n t , S n t and D ) corresponding to the analysis of NF for the JJI-IMZ path during the time period 20 March 2016–15 April 2016 for four threshold values (a) 4, (b) 7, (c) 11 and (d) 14, respectively. Figure format is similar to Figure 3.
Entropy 20 00199 g004
Figure 5. NT analysis of D P ; variations of the NT analysis parameters ( κ 1 , S n t , S n t and D ); JJI-KMK path, two different thresholds, (a) 0.30 and (b) 0.72, respectively; JJI-STU path, two different thresholds (c) 0.12 and (d) 0.54, respectively. Figure format is similar to Figure 3.
Figure 5. NT analysis of D P ; variations of the NT analysis parameters ( κ 1 , S n t , S n t and D ); JJI-KMK path, two different thresholds, (a) 0.30 and (b) 0.72, respectively; JJI-STU path, two different thresholds (c) 0.12 and (d) 0.54, respectively. Figure format is similar to Figure 3.
Entropy 20 00199 g005
Figure 6. NT analysis of T R ; variations of the NT analysis parameters ( κ 1 , S n t , S n t and D ); (a) an example (threshold = 2.1) of the JJI-KTU path; (b) an example (threshold = 1.5) of the JJI-STU path; JJI-AKT path, two different thresholds (c) 0.24 and (d) 0.48, respectively. Figure format is similar to Figure 3.
Figure 6. NT analysis of T R ; variations of the NT analysis parameters ( κ 1 , S n t , S n t and D ); (a) an example (threshold = 2.1) of the JJI-KTU path; (b) an example (threshold = 1.5) of the JJI-STU path; JJI-AKT path, two different thresholds (c) 0.24 and (d) 0.48, respectively. Figure format is similar to Figure 3.
Entropy 20 00199 g006
Figure 7. NT analysis examples for cases marginally satisfying criticality criteria; variations of the NT analysis parameters ( κ 1 , S n t , S n t and D ). (a) NF for the JJI-THY path (threshold = 9); (b) TR for the JJI-NSB path (threshold = 0.33); (c) DP for the JJI-AKT path (threshold = 0.10); (d) DP for the JJI-ANA path (threshold = 0.99). Figure format is similar to Figure 3.
Figure 7. NT analysis examples for cases marginally satisfying criticality criteria; variations of the NT analysis parameters ( κ 1 , S n t , S n t and D ). (a) NF for the JJI-THY path (threshold = 9); (b) TR for the JJI-NSB path (threshold = 0.33); (c) DP for the JJI-AKT path (threshold = 0.10); (d) DP for the JJI-ANA path (threshold = 0.99). Figure format is similar to Figure 3.
Entropy 20 00199 g007
Figure 8. Temporal evolutions of criticality as revealed by NT analysis for the VLF subionospheric propagation parameters N F , T R , and D P (top panel), and for the ULF magnetic field characteristics F h , F z , D h and δ D e p (bottom panel). Y-axis has no scale or units; only dates of criticality are indicated, while the differences in the y-axis dimension just serve as a way for easy visual discrimination of the criticality dates for the different parameters/characteristics, especially for overlapping/adjacent cases. The star symbols indicate the time of occurrence of the Kumamoto EQs.
Figure 8. Temporal evolutions of criticality as revealed by NT analysis for the VLF subionospheric propagation parameters N F , T R , and D P (top panel), and for the ULF magnetic field characteristics F h , F z , D h and δ D e p (bottom panel). Y-axis has no scale or units; only dates of criticality are indicated, while the differences in the y-axis dimension just serve as a way for easy visual discrimination of the criticality dates for the different parameters/characteristics, especially for overlapping/adjacent cases. The star symbols indicate the time of occurrence of the Kumamoto EQs.
Entropy 20 00199 g008
Table 1. VLF/LF transmitters and receivers information.
Table 1. VLF/LF transmitters and receivers information.
TransmitterReceiver
Call SignLocation (Frequency)Call SignLocation
JJIMiyazaki (22.2 kHz)AKTAkita
JJYFukushima (40 kHz)ANAAnan
NLKSeattle (24.8 kHz)IMZImizu
NPMHawaii (21.4 kHz)KMKKamakura
NWCW. Australia (19.8 kHz)KTUKatsuura
NSBNakashibetsu
STUSuttsu
TYHToyohashi
Table 2. NT Analysis results for the paths between all eight VLF/LF receivers and JJI transmitter. Only the receiver data which presented criticality are mentioned. Bold fonts indicate clear criticality indications, while normal fonts denote marginal indications (e.g., cases for which criticality is found for a few threshold values only, or for which the NT analysis parameters ( κ 1 , S n t , S n t and D ) marginally satisfy criticality criteria).
Table 2. NT Analysis results for the paths between all eight VLF/LF receivers and JJI transmitter. Only the receiver data which presented criticality are mentioned. Bold fonts indicate clear criticality indications, while normal fonts denote marginal indications (e.g., cases for which criticality is found for a few threshold values only, or for which the NT analysis parameters ( κ 1 , S n t , S n t and D ) marginally satisfy criticality criteria).
DateNFTRDP
28 March 2016 TYHANA
29 March 2016KMK, TYH ANA
30 March 2016 ANA
31 March 2016KTU
1 April 2016KTUKTU, KMK, TYH
2 April 2016 KTU
3 April 2016
4 April 2016
5 April 2016 KTU, AKT
6 April 2016NSBNSBKTU, ANA
7 April 2016NSB KTU, NSB
8 April 2016 STUSTU, NSB
9 April 2016KTUSTUSTU, IMZ
10 April 2016AKTSTUKMK, STU
11 April 2016IMZKTU, STU, IMZKMK, STU
12 April 2016IMZKTU, AKT, IMZKMK, AKT
13 April 2016 KTU, AKT
14 April 2016
15 April 2016KMKKMK

Share and Cite

MDPI and ACS Style

Potirakis, S.M.; Asano, T.; Hayakawa, M. Criticality Analysis of the Lower Ionosphere Perturbations Prior to the 2016 Kumamoto (Japan) Earthquakes as Based on VLF Electromagnetic Wave Propagation Data Observed at Multiple Stations. Entropy 2018, 20, 199. https://doi.org/10.3390/e20030199

AMA Style

Potirakis SM, Asano T, Hayakawa M. Criticality Analysis of the Lower Ionosphere Perturbations Prior to the 2016 Kumamoto (Japan) Earthquakes as Based on VLF Electromagnetic Wave Propagation Data Observed at Multiple Stations. Entropy. 2018; 20(3):199. https://doi.org/10.3390/e20030199

Chicago/Turabian Style

Potirakis, Stelios M., Tomokazu Asano, and Masashi Hayakawa. 2018. "Criticality Analysis of the Lower Ionosphere Perturbations Prior to the 2016 Kumamoto (Japan) Earthquakes as Based on VLF Electromagnetic Wave Propagation Data Observed at Multiple Stations" Entropy 20, no. 3: 199. https://doi.org/10.3390/e20030199

APA Style

Potirakis, S. M., Asano, T., & Hayakawa, M. (2018). Criticality Analysis of the Lower Ionosphere Perturbations Prior to the 2016 Kumamoto (Japan) Earthquakes as Based on VLF Electromagnetic Wave Propagation Data Observed at Multiple Stations. Entropy, 20(3), 199. https://doi.org/10.3390/e20030199

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop