Next Article in Journal
Quasi-Linearization Analysis for Entropy Generation in MHD Mixed-Convection Flow of Casson Nanofluid over Nonlinear Stretching Sheet with Arrhenius Activation Energy
Next Article in Special Issue
OpenFOAMTM Simulation of the Shock Wave Reflection in Unsteady Flow
Previous Article in Journal
A Text Classification Model via Multi-Level Semantic Features
Previous Article in Special Issue
Poisson Stability in Symmetrical Impulsive Shunting Inhibitory Cellular Neural Networks with Generalized Piecewise Constant Argument
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Critical Dynamics in Stratospheric Potential Energy Variations Prior to Significant (M > 6.7) Earthquakes

by
Dimitrios Z. Politis
1,
Stelios M. Potirakis
1,2,*,
Subrata Kundu
3,
Swati Chowdhury
3,
Sudipta Sasmal
4 and
Masashi Hayakawa
5,6
1
Department of Electrical and Electronics Engineering, Ancient Olive Grove Campus, University of West Attica, GR-12244 Egaleo, Greece
2
Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, National Observatory of Athens, Metaxa and Vasileos Pavlou, Penteli, GR-15236 Athens, Greece
3
Department of Ionospheric Sciences, Indian Center for Space Physics, 43, Chalantika, Garia 6 Station Road, Kolkata 700084, India
4
Institute of Astronomy, Space and Earth Science, AJ 316, Sector II, Salt Lake, Kolkata 700091, India
5
Hayakawa Institute of Seismo Electromagnetics, Co., Ltd., UEC Alliance Center, Chofu 11, Tokyo 182-0026, Japan
6
Advanced Wireless Communication Research Center, University of Electro-Communications (UEC), 1-5-1 Chofugaoka, Chofu, Tokyo 182-8585, Japan
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(9), 1939; https://doi.org/10.3390/sym14091939
Submission received: 26 August 2022 / Revised: 8 September 2022 / Accepted: 14 September 2022 / Published: 18 September 2022
(This article belongs to the Special Issue Symmetry in Nonlinear Dynamics and Chaos)

Abstract

:
Lithosphere–atmosphere–ionosphere coupling (LAIC) is studied through various physical or chemical quantities, obtained from different sources, which are observables of the involved complex processes. LAIC has been proposed to be achieved through three major channels: the chemical, the acoustic, and the electromagnetic. Accumulated evidence supporting the acoustic channel hypothesis has been published, while atmospheric gravity waves (AGWs) play a key role in LAIC as the leading mechanism for the transmission of energy from the lower atmosphere to the stratosphere and mesosphere, associated with atmospheric disturbances observed prior to strong earthquakes (EQs). The seismogenic AGW is the result of temperature disturbances, usually studied through stratospheric potential energy (EP). In this work, we examined 11 cases of significant EQs (M > 6.7) that occurred during the last 10 years at different geographic areas by analyzing the temperature profile at the wider location of each one of the examined EQs. The “Sounding of the Atmosphere using Broadband Emission Radiometry” (SABER) instrument, part of the “Thermosphere Ionosphere Mesosphere Energetics Dynamics” (TIMED) satellite, data were employed to compute the potential energy (EP) of the AGW. Using the temperature profile, we first calculated EP and determined the altitudes’ range for which prominent pre-seismic disturbances were observed. Subsequently, the EP time series at specific altitudes, within the determined “disturbed” range, were for the first time analyzed using the criticality analysis method termed the “natural time” (NT) method in order to find any evidence of an approach to a critical state (during a phase transition from a symmetric phase to a low symmetry phase) prior to the EQ occurrence. Our results show criticality indications in the fluctuation of EP a few days (1 to 15 days) prior to the examined EQs, except from one case. In our study, we also examined all of the temperature-related extreme phenomena that have occurred near the examined geographic areas, in order to take into account any possible non-seismic influence on the obtained results.

1. Introduction

There are different kinds of pre-seismic anomalies that have been observed on various physical or chemical parameters and they have been recorded from the ground up to space [1,2,3,4,5,6,7,8,9,10]. Several of these parameters have been extensively used together in a multi-parametric approach by many researchers in order to understand the mechanism(s) governing the so-called lithosphere–atmosphere–ionosphere coupling (LAIC) [11,12,13,14,15,16]. Specifically, five different hypotheses have been suggested, the so-called LAIC channels (acoustic, electromagnetic, thermal, electrostatic and chemical), providing different coupling mechanisms [17]. In this research, we focus on the acoustic channel hypothesis.
In the acoustic channel, atmospheric gravity waves (AGWs) are the key agents of LAIC and appear as atmospheric oscillations at stratospheric heights near the epicentral zone of the earthquake (EQ) during the preparation process. According to fluid dynamics, gravity waves are waves that are produced in a fluid medium or in the interface between two media, when the gravitational force or buoyancy tries to stabilize the equilibrium on them [18,19]. AGWs have the feature to propagate in an upward direction, disturbing the lower ionosphere [20]. In the middle of this propagation, specifically in the stratosphere, gravity waves ranging from the Brunt–Väisälä period to the inertial period affect the stratospheric wind, temperature, and pressure factors and are energized by convection systems, jet streams, and fronts [20]. The transported AGWs have energy and momentum that can reach the stratosphere and mesosphere, thus affecting the lower ionosphere. Fluctuations in pressure and temperature due to meteorological, seismic, geomagnetic activity, and solar terminators can cause AGWs in the atmosphere. Before EQs, AGW oscillations have been observed to appear around the EQ epicenter [21].
The existence of AGWs before EQs was first mentioned by Garmash et al. [22], Linkov et al. [23], and Shalimov [24], who also provided evidence of seismo–ionospheric coupling. Subsequently, several authors have studied AGWs before EQs by using both satellite and ground-based observations. Specifically, Miyaki et al. [25] first reported AGWs in very-low frequency (VLF) signals prior to EQs. Using data of the surface atmospheric pressure and magnetic field, Korepanov et al. [26] supported that AGW is a useful parameter for seismo–ionospheric coupling. Using a combined analysis of surface pressure, ionospheric perturbations, and ground-based ultra-low frequency (ULF) variations, Nakamura et al. [27] investigated the seismogenic effect of AGW for a few EQs by applying wavelet analysis. Their results revealed anomalies of wave-like structures over a period of 10 to 100 min for the 2004 Niigata-Chuetsu EQ by analyzing these parameters. Furthermore, Biswas et al. [28], using nighttime VLF fluctuation data and temperature profile data from the “Sounding of the Atmosphere using Broadband Emission Radiometry” (SABER) instrument, investigated the AGW activity before the Imphal EQ. Specifically, their results revealed enhancements in potential energy, E P , around the EQ epicenter one week before the occurrence of the EQ. Additionally, in another study, using the temperature profile from the ERA5 dataset and global positioning total electron content (GPS-TEC), Piersanti et al. [16] observed the presence of AGW in E P and found anomalous activity in GPS-TEC (using wavelet analysis) on the day of the 2018 Bayan EQ. Very recently, Sasmal et al. [12] examined the case of the 2020 Samos EQ using SABER temperature profile and GPS-TEC data. They observed abnormal AGW activity in E P and GPS-TEC (using wavelet analysis), 6 days and 11 days before the 2020 Samos EQ, respectively.
In the present work, we used the temperature profile data from the SABER instrument, which is a part of the “Thermosphere Ionosphere Mesosphere Energetics Dynamics” (TIMED) satellite, in order to study 11 significant EQs that occurred in different geographic areas (two of them occurred in the same area in very close dates). Kundu et al. [29] examined seven of these EQs (see Table 1) by using the SABER temperature profile and observed enhancements in E P that were attributed to AGW activity and appeared a few days before each EQ. Moreover, they investigated other atmospheric phenomena such as thunderstorms, cyclones, and effects from Kelvin waves (KWs) in the equatorial region around the time of occurrence of each EQ, in order to check for any possible contamination effect to their results due to non-EQ phenomena.
In the present paper, we investigated the potential energy data associated with these 11 EQs (two of them were considered together due to spatial and temporal vicinity) for possible indications of critical state by means of the so-called natural time (NT) analysis method [10]. NT analysis is a well-established criticality analysis method that has successfully been applied to foreshock seismicity and seismo-electromagnetic signals such as seismic electric signals (SES), ultra-low frequency (ULF) magnetic field variations, MHz fracto-electromagnetic emissions, and VLF sub-ionospheric propagation variations [6,10,30,31,32,33,34,35,36]. Recently, NT analysis has been used for the thermal channel on surface latent heat flux (SLHF) data, showing criticality indications a few days prior to the 2015 Nepal EQs and the 2016 Kumamoto EQ [37]. Since criticality has already been detected in a number of observables of EQ-related processes, it might be reasonable to investigate whether the stratospheric potential energy is also approaching a critical state prior to significant EQs. In this study, NT analysis was, for the first time, applied to an appropriately constructed E P time series, revealing criticality indications 1 to 15 days prior to all of the examined EQs except from one case. Additionally, for each studied EQ, we investigated whether any extreme atmospheric phenomena occurred over the analyzed geographic area after the date that E P was found to approach criticality to check whether the identified criticality could possibly be attributed to the preparation of extreme atmosphere-influencing phenomena other than EQs.
The rest of this article is summarized as follows. In Section 2, we provide information about the studied EQs. Section 3 comprises two subsections, the first presenting a brief overview of the computation of E P based on the temperature profile, which was obtained from SABER instrument of TIMED satellite, while the second presents the key notions for NT analysis method as applied to the E P time series. In Section 4, we present the analysis results, whereas in Section 5, we summarize the conclusions.

2. Studied EQs

In this article, we studied 11 significant EQs. Table 1 contains all of the information about these main events. It is noted that in this study, we focused only on the mainshocks and not any strong aftershocks. The earthquake preparation zone (EPZ) radius for each of these EQs was calculated by using the Dobrovolsky’s radius formula as: = 10 0.43 M , where M is the magnitude of the EQ [38], whereas the critical zone (CZ) radius of each EQ was calculated as: R c r = 10 0.44 M 0.78 [39], both included in Table 1. Hereafter, each EQ will be referred to using the corresponding epicenter location name, as appears in Table 1. Figure 1 shows a world map presenting all of the studied EQs epicenters as well as their EPZ and CZ.

3. Methods and Data

In this section, we briefly present the preprocessing of the SABER temperature profile data for the construction of a potential energy ( E P ) time series that can be analyzed using the NT analysis as well as key notions of the NT analysis method. Specifically, in Section 3.1 we present the computations of E P and the conversion into one time series per altitude for a specific studied geographic area, whereas in Section 3.2, we present key information about the NT analysis method as well as the basic steps of its application to the resulting E P time series.

3.1. Computation of Potential Energy (EP) Using SABER/TIMED Temperature Profile

The TIMED satellite was launched by the National Aeronautics and Space Agency (NASA) on 7 December 2001 in order to study the physical and chemical processes acting within and upon the coupled mesosphere, lower thermosphere, and ionosphere at altitudes from 60 to 180 km. The satellite is at 625 km orbit with a 74.1° inclination and includes four main instruments [40]. These are the Global Ultraviolet Imager (GUVI), the Solar Extreme Ultraviolet Experiment (SEE), the TIMED Doppler Interferometer (TIDI), and the SABER. In our investigation, we used SABER data. The specific instrument measures the ambient temperature profile for the altitudes ranging from 20 to 100 km and has a wavelength range from 1.27 to 17 μm. The geographic region that the SABER instrument can scan lies within the latitudinal range of 50° S–82° N from north viewing, while from the southward view is 50° N–82° S. Additionally, it updates its coverage every 60 days by changing its observation viewing technique. The procedure to obtain the temperature profile from the SABER instrument recordings can be found in [41].
Many researchers have used the temperature profile in order to examine the AGW activity prior to EQs by following different methods (e.g., [12,16,42,43,44,45]). In this study, we followed the procedure described in [12,28,29]. First, we took the temperature profile data from the SABER archive (http://saber.gats-inc.com/) (accessed on 1 August of 2022) in the altitude range of 20–100 km for the region of interest that includes the EPZ of the examined EQ. Subsequently, we calculated the logarithm of each individual temperature profile. After that, a third-degree polynomial was fitted to the logarithmic temperature profile. Then, all the residual values were computed by subtracting the fitted temperature profile from the initial profile. For the reason that the AGW have wavelengths longer than 4 km, a boxcar filter was applied to the residuals of the individual profiles in order to remove other wavelengths that were shorter than 4 km. The filtered residual values were set back together with the fitted profile and provide the final profile. The antilogarithm of the final fitted profile was next calculated to obtain the daily zonal mean temperature and other zonal wave components from 1 to 5. Thus, the acquired background temperature profile ( T 0 ) was computed by the summation of all of the wave components from 0 to 5. Subsequently, the perturbation temperature ( T P ) was calculated by the subtraction of the background temperature profile ( T 0 ) from the initial profiles. Finally, the potential energy, E P , which is related to AGW, is calculated using the obtained values of T P and T 0 as:
E P = 1 2 g N 2 T P T 0 2 ,
where g denotes the gravity acceleration and N is the Brunt–Väisälä frequency defined as:
N 2 = g T 0 T 0 z + g c p   ,
where z is the altitude and cp is the specific heat at a constant pressure.
Upon the completion of the calculation of E P , a nine-dimensional matrix was obtained containing all of the information involved in the calculation. Specifically, the geographic coordinates (latitude and longitude), date or day of the year (DOY), time, original SABER temperature profile, reconstructed temperature profile, perturbation temperature, Brunt–Väisälä frequency, and AGW related potential energy ( E P ) were included.
In order to apply the NT analysis to the E P data, appropriate one-dimensional signals have to be constructed. Since this was the first time NT analysis was applied to E P data, we decided to focus our investigation on the time that criticality is reached and the altitude at which the E P reaches criticality. Therefore, one time series per altitude, containing the daily E P values (one E P value per day) for each studied geographic area was constructed. The per-altitude daily-valued time series of E P were built by the following procedure.
First, we re-organized the E P data into separate datasets, one per day, containing the spatial (latitude, longitude, and altitude) distribution of E P . Since significant E P anomalies have already been identified at specific altitudes or within specific altitude ranges (third column of Table 2), for each specific EQ, we focused on these altitudes. For example, since for the Tohoku EQ the maximum E P appeared at 43 km, our analysis focused on altitudes in the range of 42.04 km–44.05 km. For each day, we re-sampled/interpolated in the altitude scale at a 0.067 km resolution by calculating the mean value of all of the E P data corresponding to the whole geographic area of interest that fall within each 0.067 km altitude zone. Finally, we constructed the final daily time series of E P (value of E P per day) for each altitude of interest (e.g., in the case of Tohoku EQ for 42.040 km, 42.107 km, …, 44.050 km). Table 2 presents the information about the enhancement of E P before the examined EQs as obtained from different already published articles except for the Miyako and Fukushima EQs, for which the results are under review.
Note that due to their spatial and (very close) temporal vicinity, the Indian Ocean−1 and Indian Ocean−2 EQs were considered as one case in our analysis.

3.2. Natural Time (NT) Analysis Method

The NT analysis method was first proposed for the analysis of ultra-low frequency (≤1 Hz) seismic electric signals (SES) [46,47,48] and has been shown to be optimal for enhancing the signals in the time–frequency space [49]. NT is a time series analysis method that can reveal when a complex system approaches critical state [10]. The wide variety of applications of NT include the analysis of foreshock seismicity, seismo-electromagnetic signals, and electrocardiograms that have extensively been explored during the last two decades [6,10,30,31,32,50,51,52]. Other newly appeared applications of NT analysis refer to the elucidation of LAIC by the analysis of various observables [37,53,54]. In this section, we present a brief overview of the key notions of the NT method, as well as its application to the E P time series. Readers are referred to [10] for a full theoretical presentation of NT analysis, and to [6] for a detail description of the application of NT to foreshock seismicity and various seismo-electromagnetic signals.
Taking a time series of L events and ignoring their “timestamp”, we define the k-th event’s “natural time” (NT),   χ k , as:
χ k = k / L
Simply put, the NT of an event, in a time series of L events, is its order of occurrence normalized to the interval (0,1] (see Figure 2). The next step is to determine the “energy” of each event in NT, which is denoted as Q k for the k-th event. It is noted that the “energy”, Q k , corresponds to different physical quantities (and not necessarily to energy) for different time series [10]. For example, the Q k of each event in the case of seismicity is the energy release (seismic moment), while for dichotomous SES signals, Q k corresponds to the SES pulse duration [10,47], etc. However, in the case of the MHz-kHz fracto-electromagnetic emissions, which are non-dichotomous signals, Q k corresponds to the energy of each event by using consecutive amplitude values above a noise threshold, as described in [32].
Subsequently, we study the new time series χ k ,   Q k Specifically, in order to identify the approach of the dynamical system to criticality, we pay attention to the following quantities: the variance, κ 1 ,
κ 1 = k = 1 L p k χ k 2 k = 1 L p k χ k 2 ,
where p k = Q k n = 1 L Q n is the normalized energy released during the k-th event, and the entropy in NT, S n t ,
S n t χ ln χ χ ln χ = k = 1 L p k χ k ln χ k k = 1 L p k χ k ln k = 1 L p k χ k ,
by considering f χ = k = 1 L p k f χ k . The entropy in NT, S n t , depends on the order of the events (i.e., it is a dynamic entropy [10,55]). The entropy in NT under time reversal, S n t , is calculated in the same manner as S n t (i.e., using Equation (5), but by analyzing in ΝΤ, the time series obtained upon considering the time reversal of the original time series). Specifically, if T ^ is the time reversal operator, its effect on a series of L events Q k ,   k = 1 , 2 , , L , is such that:
T ^ Q k = Q L k + 1 ,
which means that the first event ( k = 1 ) is positioned last in the time reversed time-series, the second becomes last but one, etc. (see Figure 2) [10,55]. Thus, the time reversal operator T ^ acting on the normalized energy, p k , results in:
T ^ p k = p L k + 1 ,
and thus S n t is, in general, different from S n t , since the natural time χ k in each time series case, represents the considered order of events (see Figure 2) [10,55].
In many dynamical complex systems, it has been found that the value of κ 1 is a measure to quantify the extent of the organization of the system at the onset of critical state [10]. The criticality is reached when (i) κ 1 has the value κ 1 = 0.07 , and (ii) 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 [10,56], where S u is the entropy of a “uniform” distribution in NT [10,55].
In the special case of the NT analysis of foreshock seismicity [10,47,48,55], we study the evolution of the quantities κ 1 , S n t , S n t , and D with time, where D is the “average” distance between the normalized power spectrum ω ˜ = k = 1 L p k e x p j ω ˜ χ k 2 , ( ω ˜ stands for the “angular frequency in NT”) of the evolving seismicity, and the theoretical estimation of ω ˜ for κ 1 = 0.07 , c r i t i c a l ω ˜ 1 κ 1 ω ˜ 2 . Moreover, an “event” for the NT analysis is considered to be any time series value of the original seismicity time series that surpasses a magnitude threshold, M t h r e s .
The analysis starts with an appropriate low threshold, and by taking into account only an adequate number of the events, which are first in the order of the occurrence. Next, the subsequent events, in the original order, are one-by-one taken into account. For each additional event that is taken into account, the quantity χ k is rescaled within the interval (0,1], while the normalized energy p k and the values of κ 1 , S n t , S n t , and D are all re-calculated. In this way, the temporal evolution of these quantities is obtained by taking into account all of the preceding events and the current event each time, until the last event has been included in the analysis. The described procedure is repeated for several increasing values of M t h r e s for each studied geographic area and everything is repeated for different overlapping areas.
The foreshock seismicity is considered to be in a true critical state, or as it is often mentioned, a “true coincidence” is achieved as soon as (i) κ 1 takes the value κ 1 = 0.07 ; (ii) the two entropies S n t , S n t satisfy the condition S n t , S n t < S u , while at the same time the following three additional conditions are satisfied: (iii) the “average” distance D should be smaller than 10 2 (i.e., D = ω ˜ c r i t i c a l ω ˜ < 10 2 , which is a practical criterion that signals the achievement of spectral coincidence) [10]; (iv) the parameter κ 1 should approach the value κ 1 = 0.07 “by descending from above” (i.e., before the main event, the parameter should gradually decrease until it reaches the critical value of 0.07; this rule was found empirically) [10,48]); and (v) the above-mentioned conditions (i)–(iv) should continue to be satisfied, even if the considered M t h r e s or the area within which the seismicity is studied are changed (within reasonable limits).
The use of the magnitude threshold excludes some of the weaker EQ events (those events whose magnitude is < M t h r e s ) from the NT analysis. However, the usage of the magnitude threshold is valid for the reason that for each seismographic network, a specific magnitude threshold, M c , is usually defined to assure data completeness, while magnitudes below M c are not considered reliable. On the other hand, the application of various M t h r e s values is useful in determining the time range within which criticality is reached. This is because, in some cases, it can be found that multiple time points may satisfy the rest of the NT critical state conditions (i)–(iv), and criterion (v) is the one that finally reveals the true time of criticality.
For the analysis of the constructed, per-altitude, daily-valued, E P time series, we applied the NT analysis following the same procedure used for the analysis of seismicity, as described above. As the “energy”, Q k , of the k-th event, we directly considered the per-altitude E P time series values without any kind of preprocessing, since the analyzed daily E P values correspond to energy. Then, the four parameters ( κ 1 ,   S n t ,   S n t ,   D ) were progressively calculated, as described above, by defining a variable threshold E P _ t h r e s and taking into account for the calculations only the daily values E P > E P _ t h r e s . Subsequently, we repeated the same procedure many times, calculating simultaneously each time the above-mentioned parameters by increasing the threshold.

4. Analysis Results

As already mentioned in the introduction, the application of NT analysis to the data of stratospheric potential energy ( E P ) is presented here for the first time. As explained in Section 3.2, the NT analysis was applied to per-altitude, daily-valued, E P time series (see Section 3.1), similar to the way it is applied to foreshock seismicity. In this work, we analyzed the E P data at different altitudes prior to 11 strong EQs that have occurred worldwide (see Section 2), in order to check whether there is at least one altitude for which the criticality was approached a few days prior to each mainshock. For each case, a long enough period of time (>30 days) before the EQ including the date of occurrence of the EQ was analyzed.
At this point, we have to clarify that for each EQ, we carefully checked whether any extreme atmospheric phenomena occurred over the analyzed geographic area after the date that E P was found to approach criticality. In the case that indications of criticality are identified before an EQ, while, a few days after the criticality is reached, an extreme atmospheric phenomenon that could affect E P also happened over the examined geographic area, this means that the criticality could well be attributed to either the pre-seismic LAIC procedure or to the process organizing the extreme atmospheric phenomenon. Therefore, in such a case, the identified criticality cannot be definitively characterized as EQ-related, so in the following, we characterize such results as “contaminated” by other atmosphere-influencing phenomena.
We checked for typhoons, convective system formations, and thunderstorms as well as for significant increases in atmospheric pressure. All of the information on the meteorological conditions were gathered from the Japan Meteorological Agency (JMA) (http://www.jma.go.jp/jma/index.html) (accessed on 1 August of 2022), the Bureau of Meteorology (Australian Government) (http://www.bom.gov.au) (accessed on 1 August 2022), the Indian Meteorological Department (IMD) (https://metnet.imd.gov.in/imdnews/) (accessed on 1 August 2022), and the joint collaboration of National Weather Service (NWS) and National Oceanic and Atmospheric Administration (NOAA) (https://www.nhc.noaa.gov/data/tcr/) (accessed on 1 August 2022). Furthermore, we also marked equatorial KWs activity, which play a vital role in AGW excitation and contaminate the E P time series for some of the examined EQs that are located near the equinox. Finally, it is noted that information about the meteorological phenomena and effects from equatorial KWs for some of the studied EQs has also been obtained from Kundu et al. [29].
Taking into account all of the above-mentioned aspects regarding the NT analysis of the per-altitude, daily-valued, E P time series as well as any other important information that is related with the analysis, we first present an example of NT analysis for the Fukushima EQ (cf. Table 1). Figure 3 presents the NT analysis results for the temporal variation of the daily E P values at the altitude of 43.516 km during the time period from 21 October 2016 until 01 December 2016 (see Figure 4). As one can see from Figure 3, all criticality conditions were satisfied for the date 14 September 2016 (i.e., 7 days prior to the Fukushima EQ (21 September 2016)). The critical state is achieved by the system due to the AGW activity (expressing the acoustic channel of LAIC) during the EQ preparation, while, as expected, the system departs from the critical state before the EQ occurrence. The NT analysis of the E P time series at different altitudes showed that criticality was approached at several altitudes in dates ranging from 6 to 9 days before the EQ occurrence (see Table 3). Moreover, we should mention that no meteorological phenomenon that could be related to the revealed criticality took place during that time period, meaning that this approach to criticality can only be associated with the pre-seismic LAIC phenomenon.
Subsequently, we present another interesting example of NT analysis of the per-altitude daily E P concerning two EQs where their epicenters were relatively close (taking into account the spatial accuracy of the used satellite data), so the same geographic region was considered for both of them, while their occurrence times were 17 days apart. These EQs were the two Nepal EQs: Ghorka−1 and Ghorka−2 (cf. Table 1).
By observing Figure 5, it is profound that the criticality criteria were clearly satisfied on 18–19 April 2015, ~1 week prior to the Ghorka−1 EQ (25 April 2015) for the 37.087 km altitude. As expected, the system departs from critical state before the Ghorka−1 EQ occurrence. However, at the specific altitude, there was no following indication of criticality that could be attributed to the Ghorka−2 EQ (12 May 2015). This is probably due to the fact that the fluctuation of E P at 37.087 km is dominated by information related to Ghorka−1 EQ, so that the embedded critical dynamics related to the Ghorka−1 EQ are “masking” any critical dynamics in the E P fluctuation that might be related to the following Ghorka−2 EQ. Therefore, the only way to find any indication of critical dynamics possibly related to the Ghorka−2 EQ is to examine other altitudes by analyzing the corresponding potential energy time series.
Indeed, Figure 6 shows that at a significantly lower altitude (at 33.010 km), the NT analysis criticality criteria were satisfied from 17 April 2015 until 1 May 2015, which is an unusually prolonged period of criticality (15 days), starting before the Ghorka−1 EQ (actually containing the time period that criticality was found for the altitude of 37.087 km) and ending 11 days prior to the Ghorka−2 EQ. This is an interesting result, implying that the fluctuation of E P at 33.010 km was such that it was able to provide information (imprints of critical dynamics) associated with both EQs. Specifically, the E P at 33.010 km first reached criticality due to the AGW activity (expressing the acoustic channel of LAIC) related to the preparation of the Ghorka−1 EQ, but the system did not depart from the critical state before the Ghorka−1 EQ occurrence because, at that altitude, the AGW activity related to the preparation of the Ghorka−2 EQ kept the system at critical state until 11 days before the Ghorka−2 EQ. Continuing our investigation by analyzing the E P time series at more altitudes within the range of interest, it was concluded that criticality indications that could be attributed to the preparation of the Ghorka−2 EQ were found at different altitudes for dates ranging from 26 April 2015 to 04 May 2015 (i.e., up to 8 days before the EQ, see Table 3).
Finally, it is noted that there were no significant meteorological phenomena over the studied area that could influence the dynamics of E P during the studied period. For this reason, the above-mentioned results were considered as definitely EQ related.
Another interesting example of NT analysis for the per-altitude daily-valued potential energy time series concerns the Chiapas EQ (cf. Table 1). The NT analysis results of E P at 45.995 km from 08 August 2017 until 18 September 2017 are presented in Figure 7. From Figure 7, it is evident that criticality was clearly reached on 29 August 2017, 10 days prior to the Chiapas EQ (08 September 2017). However, it should be mentioned that a hurricane, “Katia”, passed over the studied area during the examined time period (cf. Table 3), specifically from 05 September 2017 to 09 September 2017, a few days after the system approached criticality. Therefore, the fluctuations of E P in the examined geographic area during the days preceding the hurricane might also be influenced by the preparation of this extreme atmospheric event. Consequently, the identified criticality cannot be definitely attributed to AGW activity developed as part of LAIC (acoustic channel) during the Chiapas EQ preparation.
All of the results obtained by the NT analysis of the per-altitude daily valued E P time series for the identification of critical dynamics due to AGW activity, which are possibly related to the preparation of the studied EQs, as key agents of LAIC (acoustic channel), are summarized in Table 3. Moreover, where applicable, information about other (non-EQ-related) phenomena that might have influenced the criticality analysis results is also provided.
As one can see from Table 3, criticality was approached for all but one of the examined mainshocks. Specifically, for the Samos EQ (30 October 2020), it was not possible to find a critical signature in the E P fluctuations at any of the examined altitudes (from ~45 to ~49.5 km). At this point, it should be reminded (see also Section 3.1) that the NT analysis was applied to altitude ranges where enhancements of stratospheric potential energy, associated with AGW activity, had been identified before each of the examined EQs, as reported by different already published articles (see Table 2). However, one cannot exclude the possibility of the presence of criticality in E P at other altitudes where the organization of AGWs deploys.
Comparing the results of increased stratospheric potential energy related to AGW excitation, as presented in Table 2, with the results of the NT analysis obtained in this work (Table 3), it can be said that the time frame of appearance of the maximization of E P and the approach to criticality (last day of criticality) are both within 1–2 weeks before the EQ.
For all of the EQs that E P approached criticality, this was identified for more than one altitude, either for a continuous altitude range or for individual altitudes. Moreover, the criticality dates were not always the same at all “critical altitudes”. For the cases that the criticality date(s) changed with altitude, no specific pattern of evolution of the criticality date over altitude was identified. In other words, the criticality date does not approach the EQ date by an increase or decrease in altitude.

5. Conclusions

AGW is a well-proven phenomenon, widely studied for the elucidation of LAIC, and it is characterized as a dominant factor in the acoustic channel [21]. In this research, we applied, for the first time, the NT analysis method in the acoustic channel of LAIC in order to investigate any critical state indications prior to each of the 11 studied significant main seismic events that have happened worldwide. Specifically, the NT analysis was applied to time series of per-altitude daily values of stratospheric potential energy ( E P ), which were obtained from the temperature profiles recorded by the SABER instrument of the TIMED satellite in an altitude range of 30 to 50 km, focusing on altitudes that prominent pre-seismic disturbances were observed. In general, one cannot exclude the possibility that criticality indications can be found at various altitudes as well as different times per altitude, since AGW disturbances are upward traveling. Therefore, for each EQ case, we analyzed the per-altitude E P time series, focusing on the altitude ranges where the AGW related E P had already been identified, in order to check whether there is at least one altitude for which the criticality was approached a few days prior to the mainshock. Of course, we checked for any other (non-EQ-related) atmospheric phenomena that could possibly affect the fluctuation of E P and thus “contaminate” the obtained results.
Our results indicate that criticality was reached for all but one of the examined EQs within a period of 1–2 weeks prior to each one of them, which is consistent with various works studying the LAIC related to the specific EQs [6,30,33,34,36,57,58,59]. However, in some of the studied EQ cases, non-EQ-related phenomena that could “contaminate” the obtained results indeed existed, so for these cases, the identified criticality could not be definitively attributed to the EQ preparation processes.
This is the first work to focus on the appearance of critical signatures in one channel of LAIC (acoustic channel) for multiple EQs. However, since different channels for the LAIC phenomenon have been suggested, an interesting future research direction would be to investigate the observables of other channels for the same set of EQs in the context of the organization of the system to critical state. It is now widely accepted that understanding the LAIC mechanism calls for a simultaneous study of multi-parameters associated with the different channels of LAIC. As the seismogenic response at the lower, middle, and upper atmosphere are completely different, different sets of parameters are required to understand this subject. More specifically, we should analyze the observables of other channels, for example, for the thermal channel (surface latent heat flux, relative humidity, etc.), for the electromagnetic channel (VLF radio anomalies, total electron content (TEC), energetic particle precipitation, etc.), etc., in order to better understand which mechanism(s) is/are involved in the LAIC for each EQ case. Although the stratospheric AGW is formed due to the seismogenic perturbation of temperature, originating from a thermal modulation, the nature of the thermal excitation in the troposphere and stratosphere for a single EQ will be completely different. The original temperature, latent heat, relative humidity, etc., which is the prime source of thermal convection to build up such waves, has a different perturbation profile in the lower troposphere or at the near-Earth surface shows different seismogenic modulation. For some of the examined EQs, different observables have been analyzed in terms of criticality. For example, the depression of the horizontal magnetic field component, which is considered to be associated with ionospheric anomalies, has been analyzed for the Tohoku EQ [36] and the Kumamoto EQs [6,30], VLF sub-ionospheric propagation variations have been analyzed for the Kumamoto EQ [6,59], whereas, for the same EQ, GNSS (Global Navigation Satellite System) deformation data and its fluctuations in two AGW bands have also been analyzed [54], and very recently the SLHF (surface latent heat flux) time series for the 2016 Kumamoto EQ and 2015 Nepal EQs have also been analyzed [37]. However, a systematic investigation of the LAIC mechanism(s) calls for the analysis of different channel observables for the same set of EQs.

Author Contributions

Conceptualization, S.M.P., S.S. and M.H.; Methodology, D.Z.P., S.M.P., S.K., S.C., S.S. and M.H.; Software, S.M.P., S.K., S.C. and S.S.; Validation, D.Z.P., S.M.P., S.K., S.C., S.S. and M.H.; Formal analysis, D.Z.P., S.K. and S.C.; Investigation, D.Z.P., S.M.P., S.K. and S.C.; Data curation, D.Z.P. and S.K.; Writing—original draft preparation, D.Z.P., S.M.P. and S.S.; Writing—review and editing, D.Z.P., S.M.P., S.K., S.C., S.S. and M.H.; Visualization, D.Z.P.; Supervision, S.M.P. All authors have read and agreed to the published version of the manuscript.

Funding

The research received no external funding.

Data Availability Statement

The temperature profile data were taken from the SABER archive (http://saber.gats-inc.com/ (accessed on 1 August 2022)). The data of the constructed per-altitude time series of potential energy are available from S.K. or D.Z.P. upon request.

Acknowledgments

We are thankful to the SABER and their research staff for providing the atmospheric temperature profile. We also thank the Japan Meteorological Agency (JMA), the Japan Meteorological Agency (JMA) (http://www.jma.go.jp/jma/index.html (accessed on 1 August 2022)), the Bureau of Meteorology (Australian Government) (http://www.bom.gov.au (accessed on 1 August 2022)), the Indian Meteorological Department (IMD) (https://metnet.imd.gov.in/imdnews/ (accessed on 1 August 2022)), the joint collaboration of the National Weather Service (NWS) and the National Oceanic and Atmospheric Administration (NOAA) (https://www.nhc.noaa.gov/data/tcr/ (accessed on 1 August 2022)) for the access to their databases.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hayakawa, M.; Molchanov, O.A. Seismo Electromagnetics: Lithosphere-Atmosphere-Ionosphere Coupling; TERRAPUB: Tokyo, Japan, 2002. [Google Scholar]
  2. Molchanov, O.A.; Hayakawa, M. Seismo Electromagnetics and Related Phenomena: History and Latest Results; TERRAPUB: Tokyo, Japan, 2008. [Google Scholar]
  3. 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]
  4. Hayakawa, M. Earthquake Prediction with Radio Techniques; Wiley: Singapore, 2015. [Google Scholar]
  5. Eftaxias, K.; Potirakis, S.M.; Contoyiannis, Y. Four-Stage model of earthquake generation in terms of fracture-induced electromagnetic emissions a review. In Complexity of Seismic Time Series: Measurement and Application; Chelidze, T., Vallianatos, F., Telesca, L., Eds.; Elsevier: Oxford, UK, 2018; pp. 437–502. [Google Scholar]
  6. Potirakis, S.M.; Contoyiannis, Y.; Schekotov, A.; Eftaxias, K.; Hayakawa, M. Evidence of critical dynamics in various electromagnetic precursors. Eur. Phys. J. Spec. Top. 2021, 230, 151–177. [Google Scholar] [CrossRef]
  7. Uyeda, S.; Hayakawa, M.; Nagao, T.; Molchanov, O.; Hattori, K.; Orihara, Y.; Gotoh, K.; Akinaga, Y.; Tanaka, H. Electric and magnetic phenomena observed before the volcano-seismic activity in 2000 in the Izu Island Region, Japan. Proc. Natl. Acad. Sci. USA 2002, 99, 7352–7355. [Google Scholar] [CrossRef] [PubMed]
  8. Uyeda, S.; Nagao, T.; Kamogawa, M. Short-term earthquake prediction: Current status of seismo-electromagnetics. Tectonophysics 2009, 470, 205–213. [Google Scholar] [CrossRef]
  9. Uyeda, S. On earthquake prediction in Japan. Proc. Jpn. Acad. Ser. B 2013, 89, 391–400. [Google Scholar] [CrossRef] [PubMed]
  10. Varotsos, P.; Sarlis, N.V.; Skordas, E.S. Natural Time Analysis: The New View of Time Precursory Seismic Electric Signals, Earthquakes and Other Complex Time Series; Springer: Berlin, Germany, 2011. [Google Scholar]
  11. Hayakawa, M.; Izutsu, J.; Schekotov, A.; Yang, S.-S.; Solovieva, M.; Budilova, E. Lithosphere–atmosphere–ionosphere coupling effects based on multiparameter precursor observations for February–March 2021 earthquakes (M~7) in the offshore of Tohoku area of Japan. Geosciences 2021, 11, 481. [Google Scholar] [CrossRef]
  12. Sasmal, S.; Chowdhury, S.; Kundu, S.; Politis, D.Z.; Potirakis, S.M.; Balasis, G.; Hayakawa, M.; Chakrabarti, S.K. Pre-seismic irregularities during the 2020 Samos (Greece) earthquake (M = 6.9) as investigated from multi-parameter approach by ground and space-based techniques. Atmosphere 2021, 12, 1059. [Google Scholar] [CrossRef]
  13. Pulinets, S.; Ouzounov, D.; Davydenko, D.; Petrukhin, A. Multiparameter monitoring of short-term earthquake precursors and its physical basis. implementation in the Kamchatka Region. E3S Web Conf. 2016, 11, 00019. [Google Scholar] [CrossRef]
  14. De Santis, A.; Cianchini, G.; Marchetti, D.; Piscini, A.; Sabbagh, D.; Perrone, L.; Campuzano, S.A.; Inan, S. A multiparametric approach to study the preparation phase of the 2019 M7.1 Ridgecrest (California, United States) earthquake. Front. Earth Sci. 2020, 8, 540398. [Google Scholar] [CrossRef]
  15. Chetia, T.; Sharma, G.; Dey, C.; Raju, P.L. Multi-parametric approach for Earthquake Precursor Detection in assam valley (eastern Himalaya, India) using satellite and Ground Observation Data. Geotectonics 2020, 54, 83–96. [Google Scholar] [CrossRef]
  16. Piersanti, M.; Materassi, M.; Battiston, R.; Carbone, V.; Cicone, A.; D’Angelo, G.; Diego, P.; Ubertini, P. Magnetospheric–ionospheric–lithospheric coupling model. 1: Observations during the 5 August 2018 Bayan earthquake. Remote Sens. 2020, 12, 3299. [Google Scholar] [CrossRef]
  17. Hayakawa, M.; Schekotov, A.; Izutsu, J.; Nickolaenko, A.P. Seismogenic effects in ULF/ELF/VLF electromagnetic waves. Int. J. Electron. Appl. Res. 2019, 6, 1–86. [Google Scholar] [CrossRef]
  18. Phillips, O.M. On the generation of waves by turbulent wind. J. Fluid Mech. 1957, 2, 417. [Google Scholar] [CrossRef]
  19. Lighthill, J. Waves in Fluids; Cambridge University Press: Cambridge, UK, 2003. [Google Scholar]
  20. Fritts, D.C.; Alexander, M.J. Gravity wave dynamics and effects in the Middle Atmosphere. Rev. Geophys. 2003, 41. [Google Scholar] [CrossRef]
  21. Hayakawa, M.; Kasahara, Y.; Nakamura, T.; Hobara, Y.; Rozhnoi, A.; Solovieva, M.; Molchanov, O.; Korepanov, V. Atmospheric gravity waves as a possible candidate for seismo-ionospheric perturbations. J. Atmos. Electr. 2011, 31, 129–140. [Google Scholar] [CrossRef]
  22. Garmash, S.V.; Linkov, E.N.; Petrova, L.N.; Shved, G.N. Excitation of atmospheric oscillations by seismogravitational vibrations of the earth. Izv. Akad. Nauk. SSSR Fiz. Atmos. Okeana 1989, 35, 1290–1299. [Google Scholar]
  23. Linkov, E.M.; Petrova, L.N.; Zuroshvili, D.D. Seismogravitational vibrations of the earth and related disturbances of the atmosphere. Dokl. Akad. Nauk. SSSR 1989, 306, 315–317. [Google Scholar]
  24. Shalimov, S.L. Lithosphere-ionosphere relationship: A new way to predict earthquakes? Epis. Int. Geophys. Newsmag. 1992, 15, 252–254. [Google Scholar] [CrossRef]
  25. Miyaki, K.; Hayakawa, M.; Molchanov, O.A. The Role of Gravity Waves in the Lithosphere-Ionosphere Coupling, as Revealed from the Subionospheric LF Propagation Data. In Seismo-Electromagnetics: Lithosphere-Atmosphere-Ionosphere Coupling; Hayakawa, M., Molchanov, O., Eds.; TERRAPUB: Tokyo, Japan, 2002; pp. 229–232. [Google Scholar]
  26. Korepanov, V.; Hayakawa, M.; Yampolski, Y.; Lizunov, G. AGW as a seismo-ionospheric coupling responsible agent. Phys. Chem. Earth A/B/C 2009, 34, 485–495. [Google Scholar] [CrossRef]
  27. Nakamura, T.; Korepanov, V.; Kasahara, Y.; Hobara, Y.; Hayakawa, M. An evidence on the lithosphere-ionosphere coupling in terms of atmospheric gravity waves on the basis of a combined analysis of surface pressure, ionospheric perturbations and ground-based ULF Variations. J. Atmos. Electr. 2013, 33, 53–68. [Google Scholar] [CrossRef]
  28. Biswas, S.; Kundu, S.; Ghosh, S.; Chowdhury, S.; Yang, S.-S.; Hayakawa, M.; Chakraborty, S.K.; Chakrabarti, S.; Sasmal, S. Contaminated effect of geomagnetic storms on pre-seismic atmospheric and ionospheric anomalies during Imphal earthquake. Open J. Earthq. Res. 2020, 9, 383–402. [Google Scholar] [CrossRef]
  29. Kundu, S.; Chowdhury, S.; Ghosh, S.; Sasmal, S.; Politis, D.Z.; Potirakis, S.M.; Yang, S.-S.; Chakrabarti, S.K.; Hayakawa, M. Seismogenic anomalies in atmospheric gravity waves as observed from SABER/TIMED satellite during large earthquakes. J. Sens. 2022, 2022, 1–23. [Google Scholar] [CrossRef]
  30. 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]
  31. Potirakis, S.M.; Schekotov, A.; Contoyiannis, Y.; Balasis, G.; Koulouras, G.; Melis, N.; Boutsi, A.; Hayakawa, M.; Eftaxias, K.; Nomicos, C. On possible electromagnetic precursors to a significant earthquake (MW = 6.3) occurred in Lesvos (Greece) on 12 June 2017. Entropy 2019, 21, 241. [Google Scholar] [CrossRef]
  32. Potirakis, S.M.; Karadimitrakis, A.; Eftaxias, K. Natural time analysis of critical phenomena: The case of pre-fracture electromagnetic emissions. Chaos 2013, 23, 023117. [Google Scholar] [CrossRef]
  33. 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]
  34. Politis, D.Z.; Potirakis, S.M.; Contoyiannis, Y.F.; Biswas, S.; Sasmal, S.; Hayakawa, M. Statistical and criticality analysis of the lower ionosphere prior to the 30 October 2020 Samos (Greece) earthquake (M6.9), based on VLF electromagnetic propagation data as recorded by a new VLF/LF receiver installed in Athens (Greece). Entropy 2021, 23, 676. [Google Scholar] [CrossRef]
  35. 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]
  36. Hayakawa, M.; Schekotov, A.; Potirakis, S.M.; 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]
  37. Ghosh, S.; Chowdhury, S.; Kundu, S.; Sasmal, S.; Politis, D.Z.; Potirakis, S.M.; Hayakawa, M.; Chakraborty, S.; Chakrabarti, S.K. Unusual surface latent heat flux variations and their critical dynamics revealed before strong earthquakes. Entropy 2022, 24, 23. [Google Scholar] [CrossRef]
  38. Dobrovolsky, I.P.; Zubkov, S.I.; Miachkin, V.I. Estimation of the size of earthquake preparation zones. Pure Appl. Geophys. 1979, 117, 1025–1044. [Google Scholar] [CrossRef]
  39. Bowman, D.D.; Ouillon, G.; Sammis, C.G.; Sornette, A.; Sornette, D. An observational test of the critical earthquake concept. J. Geophys. Res. Solid Earth 1998, 103, 24359–24372. [Google Scholar] [CrossRef] [Green Version]
  40. Shuai, J.; Zhang, S.D.; Huang, C.M.; Yi, F.; Huang, K.M.; Gan, Q.; Gong, Y. Climatology of global gravity wave activity and dissipation revealed by Saber/Timed Temperature Observations. China Technol. Sci. 2014, 57, 998–1009. [Google Scholar] [CrossRef]
  41. Remsberg, E.E.; Marshall, B.T.; Garcia-Comas, M.; Krueger, D.; Lingenfelser, G.S.; Martin-Torres, J.; Mlynczak, M.G.; Russell, J.M.; Smith, A.K.; Zhao, Y.; et al. Assessment of the quality of the version 1.07 temperature-versus-pressure profiles of the middle atmosphere from timed/saber. J. Geophys. Res. Atmos. 2008, 113, D17. [Google Scholar] [CrossRef]
  42. Yan, X.; Arnold, N.; Remedios, J. Global observations of gravity waves from high resolution dynamics limb sounder temperature measurements: A year long record of temperature amplitude and vertical wavelength. J. Geophys. Res. Atmos. 2010, 115, D10. [Google Scholar] [CrossRef]
  43. Yamashita, C.; England, S.L.; Immel, T.J.; Chang, L.C. Gravity wave variations during elevated stratopause events using SABER Observations. J. Geophys. Res. Atmos. 2013, 118, 5287–5303. [Google Scholar] [CrossRef]
  44. Thurairajah, B.; Bailey, S.M.; Cullens, C.Y.; Hervig, M.E.; Russell, J.M. Gravity Wave activity during recent stratospheric sudden warming events from Sofie Temperature Measurements. J. Geophys. Res. Atmos. 2014, 119, 8091–8103. [Google Scholar] [CrossRef]
  45. Yang, S.-S.; Hayakawa, M. Gravity wave activity in the stratosphere before the 2011 Tohoku earthquake as the mechanism of lithosphere-atmosphere-ionosphere coupling. Entropy 2020, 22, 110. [Google Scholar] [CrossRef]
  46. 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]
  47. 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]
  48. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S. Spatio-temporal complexity aspects on the interrelation between seismic electric signals and seismicity. Pract. Athens Acad. 2001, 76, 294–321. [Google Scholar]
  49. 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]
  50. Baldoumas, G.; Peschos, D.; Tatsis, G.; Chronopoulos, S.K.; Christofilakis, V.; Kostarakis, P.; Varotsos, P.; Sarlis, N.V.; Skordas, E.S.; Bechlioulis, A.; et al. A prototype photoplethysmography electronic device that distinguishes congestive heart failure from healthy individuals by applying natural time analysis. Electronics 2019, 8, 1288. [Google Scholar] [CrossRef] [Green Version]
  51. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S.; Nagao, T.; Kamogawa, M. The unusual case of the ultra-deep 2015 Ogasawara earthquake (Mw 7.9): Natural Time Analysis. EPL 2021, 135, 49002. [Google Scholar] [CrossRef]
  52. Sarlis, N.V.; Skordas, E.S.; Varotsos, P.A. A remarkable change of the entropy of seismicity in natural time under time reversal before the super-giant M9 Tohoku earthquake on 11 March 2011. EPL 2018, 124, 29001. [Google Scholar] [CrossRef]
  53. Hayakawa, M.; Potirakis, S.M.; Saito, Y. Possible relation of air ion density anomalies with earthquakes and the associated precursory ionospheric perturbations: An analysis in terms of criticality. Int. J. Electron. Appl. Res. 2018, 5, 56–75. [Google Scholar] [CrossRef]
  54. Yang, S.-S.; Potirakis, S.M.; Sasmal, S.; Hayakawa, M. Natural Time Analysis of global navigation satellite system surface deformation: The case of the 2016 Kumamoto earthquakes. Entropy 2020, 22, 674. [Google Scholar] [CrossRef]
  55. Varotsos, P.A.; Sarlis, N.V.; Skordas, E.S.; Tanaka, H.K.; Lazaridou, M.S. Entropy of seismic electric signals: Analysis in natural time under time reversal. Phys. Rev. E 2006, 73, 031114. [Google Scholar] [CrossRef]
  56. Sarlis, N.V.; Skordas, E.S.; Varotsos, P.A. Similarity of fluctuations in systems exhibiting self-organized criticality. EPL 2011, 96, 28006. [Google Scholar] [CrossRef]
  57. Potirakis, S.M.; Contoyiannis, Y.; Schekotov, A.; Asano, T.; Hayakawa, M. Analysis of the ultra-low frequency magnetic field fluctuations prior to the 2016 Kumamoto (Japan) earthquakes in terms of the method of critical fluctuations. Physica A Stat. Mech. Appl. 2019, 514, 563–572. [Google Scholar] [CrossRef]
  58. Contoyiannis, Y.; Potirakis, S.M.; Eftaxias, K.; Hayakawa, M.; Schekotov, A. Intermittent criticality revealed in ULF Magnetic Fields prior to the 11 march 2011 Tohoku earthquake (Mw = 9). Physica A Stat. Mech. Appl. 2016, 452, 19–28. [Google Scholar] [CrossRef]
  59. 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. [Google Scholar] [CrossRef] [Green Version]
Figure 1. A world map presenting the epicenters along with the corresponding EPZs and CZs of the 11 studied EQs. Each black point indicates the epicenter of the corresponding EQ, while the EPZ and the CZ around each epicenter are shown as dashed and solid circles, respectively. EPZs and CZs of the Tohoku, Indian Ocean−1, Indian Ocean−2, Miyako, Gorkha−1, Gorkha−2, Kumamoto, Fukushima, Imphal, Chiapas, and Samos EQs are shown in the blue, red, light green, magenta, yellow, light blue, green, brown, purple, cyan, and orange colors, respectively.
Figure 1. A world map presenting the epicenters along with the corresponding EPZs and CZs of the 11 studied EQs. Each black point indicates the epicenter of the corresponding EQ, while the EPZ and the CZ around each epicenter are shown as dashed and solid circles, respectively. EPZs and CZs of the Tohoku, Indian Ocean−1, Indian Ocean−2, Miyako, Gorkha−1, Gorkha−2, Kumamoto, Fukushima, Imphal, Chiapas, and Samos EQs are shown in the blue, red, light green, magenta, yellow, light blue, green, brown, purple, cyan, and orange colors, respectively.
Symmetry 14 01939 g001
Figure 2. The notions of NT and time reversal through a sequence of six events of different “energy”, Q , for which the original order of occurrence is indicated by the labels “1st”, “2nd”, etc. The upper panel shows the events in the conventional time domain. One can see that these are not evenly distributed in conventional time, t . The middle panel shows the same six events in the NT domain, where the events are now evenly distributed in NT, χ . The lower panel shows the same six events after the application of the time reversal operator T ^ (see Equation (6)), where one can see that their order has been reversed and now each one of them corresponds to a different NT than in the middle panel.
Figure 2. The notions of NT and time reversal through a sequence of six events of different “energy”, Q , for which the original order of occurrence is indicated by the labels “1st”, “2nd”, etc. The upper panel shows the events in the conventional time domain. One can see that these are not evenly distributed in conventional time, t . The middle panel shows the same six events in the NT domain, where the events are now evenly distributed in NT, χ . The lower panel shows the same six events after the application of the time reversal operator T ^ (see Equation (6)), where one can see that their order has been reversed and now each one of them corresponds to a different NT than in the middle panel.
Symmetry 14 01939 g002
Figure 3. The NT analysis of the E P temporal variation at 43.516 km for the Fukushima EQ; the examined time period was from 21 October 2016 until 31 September 2016. The presented temporal variations in the NT parameters correspond to the different E P _ t h r e s thresholds: (a) 0; (b) 1.4; (c) 1.6; and (d) 1.8, respectively. The limit value for the entropy S u   0.0966 is marked as a horizontal solid light green line, while the “critical” κ 1 value 0.070 is shown as a horizontal solid grey line, along with the limits of ±0.005 around 0.070, denoted as two horizontal dashed grey lines. The 10 2   D limit is shown as a horizontal brown line. The events presented in each panel depend on the corresponding threshold. It is noted that the tick values of the x-axis show the conventional time (Date) of the occurrence of each event although the x-axis is linear to the NT (i.e., events are equally spaced in the x-axis regardless of their date of occurrence, see Figure 2).
Figure 3. The NT analysis of the E P temporal variation at 43.516 km for the Fukushima EQ; the examined time period was from 21 October 2016 until 31 September 2016. The presented temporal variations in the NT parameters correspond to the different E P _ t h r e s thresholds: (a) 0; (b) 1.4; (c) 1.6; and (d) 1.8, respectively. The limit value for the entropy S u   0.0966 is marked as a horizontal solid light green line, while the “critical” κ 1 value 0.070 is shown as a horizontal solid grey line, along with the limits of ±0.005 around 0.070, denoted as two horizontal dashed grey lines. The 10 2   D limit is shown as a horizontal brown line. The events presented in each panel depend on the corresponding threshold. It is noted that the tick values of the x-axis show the conventional time (Date) of the occurrence of each event although the x-axis is linear to the NT (i.e., events are equally spaced in the x-axis regardless of their date of occurrence, see Figure 2).
Symmetry 14 01939 g003
Figure 4. Variation in the E P in the altitude of 43.516 km for Fukushima EQ. At the specific altitude, E P was maximized on 18–19 September 2016, indicating notable AGW activity prior to the Fukushima EQ. The date of the earthquake is indicated as “EQ”.
Figure 4. Variation in the E P in the altitude of 43.516 km for Fukushima EQ. At the specific altitude, E P was maximized on 18–19 September 2016, indicating notable AGW activity prior to the Fukushima EQ. The date of the earthquake is indicated as “EQ”.
Symmetry 14 01939 g004
Figure 5. The NT analysis of E P temporal variation at 37.087 km for the Gorkha−1 EQ; the examined time period was from 10 April 2015 until 20 May 2015. The presented temporal variations of the NT parameters correspond to the different E P _ t h r e s thresholds: (a) 0; (b) 1; (c) 1.4; and (d) 1.8, respectively. The figure format follows the format of Figure 3.
Figure 5. The NT analysis of E P temporal variation at 37.087 km for the Gorkha−1 EQ; the examined time period was from 10 April 2015 until 20 May 2015. The presented temporal variations of the NT parameters correspond to the different E P _ t h r e s thresholds: (a) 0; (b) 1; (c) 1.4; and (d) 1.8, respectively. The figure format follows the format of Figure 3.
Symmetry 14 01939 g005
Figure 6. The NT analysis of E P temporal variation at 33.010 km showing criticality that can be attributed to the Gorkha−2 EQ; the examined time period was from 10 April 2015 until 20 May 2015. The presented temporal variations of the NT parameters correspond to the different E P _ t h r e s thresholds: (a) 1.6; (b) 1.7; (c) 1.9; and (d) 2, respectively. The figure format follows the format of Figure 3.
Figure 6. The NT analysis of E P temporal variation at 33.010 km showing criticality that can be attributed to the Gorkha−2 EQ; the examined time period was from 10 April 2015 until 20 May 2015. The presented temporal variations of the NT parameters correspond to the different E P _ t h r e s thresholds: (a) 1.6; (b) 1.7; (c) 1.9; and (d) 2, respectively. The figure format follows the format of Figure 3.
Symmetry 14 01939 g006
Figure 7. The NT analysis of the E P temporal variation at 45.995 km for the Chiapas EQ; the examined time period was from 08 August 2017 until 18 September 2017. The presented temporal variations in the NT parameters correspond to the different E P _ t h r e s thresholds: (a) 0; (b) 2.8; (c) 3; and (d) 3.6, respectively. The figure format follows the format of Figure 3.
Figure 7. The NT analysis of the E P temporal variation at 45.995 km for the Chiapas EQ; the examined time period was from 08 August 2017 until 18 September 2017. The presented temporal variations in the NT parameters correspond to the different E P _ t h r e s thresholds: (a) 0; (b) 2.8; (c) 3; and (d) 3.6, respectively. The figure format follows the format of Figure 3.
Symmetry 14 01939 g007
Table 1. Information about the studied EQs as provided by the United States of Geological Survey (USGS) (https://earthquake.usgs.gov/ (accessed on 1 August 2022)). The corresponding EPZ and CZ radii as well as the EQ classification according to the seven (magnitude-based) classes defined by USGS are also shown. The symbols (*), (**), and (***) in the epicenter location indicate an EQ that has been studied in [12,28,29], respectively.
Table 1. Information about the studied EQs as provided by the United States of Geological Survey (USGS) (https://earthquake.usgs.gov/ (accessed on 1 August 2022)). The corresponding EPZ and CZ radii as well as the EQ classification according to the seven (magnitude-based) classes defined by USGS are also shown. The symbols (*), (**), and (***) in the epicenter location indicate an EQ that has been studied in [12,28,29], respectively.
DateTime
(UT)
Depth
(km)
Epicenter
Latitude
Epicenter
Longitude
Magnitude
(MW)
EQ ClassEPZ
Radius
(km)
CZ
Radius
(km)
Epicenter
Location
Country/
Region
11 March 201105:46:242938.297° N142.373° E9.1Great81851674.9Tohoku *Japan
11 April 201208:38:36202.327° N93.063° E8.6Great49891009.3Indian Ocean−1 *Indo-Australian Plate
11 April 201210:43:1025.10.802° N92.463° E8.2Great3357672.9Indian Ocean−2 *Indo-Australian Plate
16 February 201523:06:282339.856° N142.881° E6.7Strong760147.2MiyakoJapan
25 April 201506:11:258.228.231° N84.731° E7.8Major2259448.7Gorkha−1 *Nepal
12 May 201507:05:191527.809° N86.066° E7.3Major1377270.3Gorkha−2 *Nepal
15 April 201616:25:061032.791° N130.754° E7.0Major1023199.5Kumamoto *Japan
21 November 201620:59:49937.393° N141.387° E6.9Strong936180.3FukushimaJapan
3 January 201623:05:225524.804° N93.651° E6.7Strong760147.2Imphal **India
8 September 201704:49:1947.415.022° N93.899° W8.2Great3357672.9Chiapas *Mexico
30 October 202011:51:272137.897° N26.784° E7.0Major1023199.5Samos ***Greece
Table 2. Information about the enhancement of EP before the examined EQs and the altitude range analyzed using NT analysis. The symbols (*), (**), and (***) in the EQ name indicate that the altitude, or the range of altitudes where significant EP enhancement appears were identified in [12,28,29], respectively.
Table 2. Information about the enhancement of EP before the examined EQs and the altitude range analyzed using NT analysis. The symbols (*), (**), and (***) in the EQ name indicate that the altitude, or the range of altitudes where significant EP enhancement appears were identified in [12,28,29], respectively.
EQ Name
and Magnitude
(MW)
EP Enhancement Time
(No. of Days before the EQ)
Maximum EP Altitude/
Altitude Range
(km)
Analyzed Altitude
Range
Tohoku *
(9.1)
4–643min: 42.040 km
step: 0.067 km
max: 44.050 km
Indian Ocean−1 & Indian Ocean−2 *
(8.6 & 8.2)
2–634–37min: 32.007 km
step: 0.067 km
max: 38.037 km
Miyako
(6.7)
9–1138–42min: 38.027 km
step: 0.067 km
max: 41.779 km
Gorkha−1 *
(7.8)
13–1435–36min: 33.000 km
step: 0.067 km
max: 37.087 km
Gorkha−2 *
(7.3)
1–334–36min: 33.010 km
step: 0.067 km
max: 37.432 km
Kumamoto *
(7.0)
8–943min: 40.033 km
step: 0.067 km
max: 44.99 km
Fukushima
(6.9)
2–444–46min: 41.037 km
step: 0.067 km
max: 45.995 km
Imphal **
(6.7)
7–1044–46min: 43.000 km
step: 0.067 km
max: 47.020 km
Chiapas *
(8.2)
742–46min: 41.037 km
step: 0.067 km
max: 45.995 km
Samos ***
(7.0)
5–746–48min: 45.117 km
step: 0.067 km
max: 49.539 km
Table 3. A summary of the NT analysis results for each EQ appearing in Table 1.
Table 3. A summary of the NT analysis results for each EQ appearing in Table 1.
EQ NameCriticality Altitude
(km)
Criticality Date
(No. of Days before EQ)
Date of EQ
Occurrence
Other Extreme
Phenomena after
the Occurrence of
Criticality
Selected Time
Period for the
Analysis
Selected Spatial
Region for the
Analysis
Tohoku43.849–44.05008 March 2011
(3)
11 March 201117–23 March 2011 were disturbed due to three low pressure cyclonic systems15 February 2011–
21 March 2011
(20°, 60°) N,
(120°, 160°) E
Indian Ocean−1 & Indian Ocean−2 *32.007–34.48607 April 2012
(4)
11 April 2012A tropical low-pressure area of 19U was formed during
16–25 April 2012. Also, equatorial KWs.
13 March 2012–
21 April 2012
(–15°S, 10° N),
(70°, 110°) E
Miyako38.228, 38.429, 38.764, 39.769, 39.836, 39.903, 39.970, 40.037, 40.104, 40.171, 40.707, 41.042, 41.77907 February 2015–
13 February 2015
(3–9)
16 February 2015-17 January 2015–
26 February 2015
(25°, 50°) N,
(130°, 160°) E
Gorkha−133.000, 33.067, 33.268, 33.469, 33.938, 34.407, 35.010, 35.144, 35.613, 36.082, 36.484, 36.752, 37.08715 April 2015–24 April 2015
(1–10)
25 April 2015-10 April 2015–
20 April 2015
(20°, 40°)N,
(70°, 110°)E
Gorkha−233.010, 33.268, 33.814, 35.556, 35.958, 36.360, 36.762, 37.030, 37.43226 April 2015–4 May 2015,
(8–16)
12 May 2015-10 April 2015–
20 May 2015
(20°, 40°)N,
(70°, 110°)E
Kumamoto40.033–40.43506 April 2016–13 April 2016
(2–9)
15 April 201614 April–19 March 2016 were disturbed due to a few frontal systems15 March 2016–
25 April 2016
(15°, 40°)N,
(120°, 140°)E
Fukushima41.908, 42.109, 42.176, 42.243, 42.980, 43.449, 43.516, 43.918, 45.660, 45.928, 45.99512 November 2016–15 November 2016
(6–9)
21 November 2016-21 October 2016–
01 December 2016
(25°, 50°)N,
(125°, 150°)E
Imphal43.000, 43.603, 44.206, 45.010, 45.479, 45.948, 46.484, 46.81926 December 2015–29 December 2015
(5–8)
03 January 2016-15 December 2016–
13 January 2016
(15°, 35°)N,
(70°, 110°)E
Chiapas45.660–45.99529 August 2017
(10)
08 September 2017Hurricane “Katia” (05–09 September 2017) in the region (22°, 25°)N, (94°, 98°)W
Also, equatorial KWs
08 August 2017–
18 September 2017
(5° S, 25° N),
(70°, 110°)W
Samos--30 October 2020-30 September 2020–
30 November 2020
(30°, 50°)N,
(10°, 40°)E
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Politis, D.Z.; Potirakis, S.M.; Kundu, S.; Chowdhury, S.; Sasmal, S.; Hayakawa, M. Critical Dynamics in Stratospheric Potential Energy Variations Prior to Significant (M > 6.7) Earthquakes. Symmetry 2022, 14, 1939. https://doi.org/10.3390/sym14091939

AMA Style

Politis DZ, Potirakis SM, Kundu S, Chowdhury S, Sasmal S, Hayakawa M. Critical Dynamics in Stratospheric Potential Energy Variations Prior to Significant (M > 6.7) Earthquakes. Symmetry. 2022; 14(9):1939. https://doi.org/10.3390/sym14091939

Chicago/Turabian Style

Politis, Dimitrios Z., Stelios M. Potirakis, Subrata Kundu, Swati Chowdhury, Sudipta Sasmal, and Masashi Hayakawa. 2022. "Critical Dynamics in Stratospheric Potential Energy Variations Prior to Significant (M > 6.7) Earthquakes" Symmetry 14, no. 9: 1939. https://doi.org/10.3390/sym14091939

APA Style

Politis, D. Z., Potirakis, S. M., Kundu, S., Chowdhury, S., Sasmal, S., & Hayakawa, M. (2022). Critical Dynamics in Stratospheric Potential Energy Variations Prior to Significant (M > 6.7) Earthquakes. Symmetry, 14(9), 1939. https://doi.org/10.3390/sym14091939

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