Investigation of Pre-Earthquake Ionospheric and Atmospheric Disturbances for Three Large Earthquakes in Mexico

: The purpose of the present study is to investigate simultaneously pre-earthquake ionospheric and atmospheric disturbances by the application of different methodologies, with the ultimate aim to detect their possible link with the impending seismic event. Three large earthquakes in Mexico are selected (8.2 Mw, 7.1 Mw and 6.6 Mw during 8 and 19 September 2017 and 21 January 2016 respectively), while ionospheric variations during the entire year 2017 prior to 37 earthquakes are also examined. In particular, Total Electron Content (TEC) retrieved from Global Navigation Satellite System (GNSS) networks and Atmospheric Chemical Potential (ACP) variations extracted from an atmospheric model are analyzed by performing statistical and spectral analysis on TEC measurements with the aid of Global Ionospheric Maps (GIMs), Ionospheric Precursor Mask (IPM) methodology and time series and regional maps of ACP. It is found that both large and short scale ionospheric anomalies occurring from few hours to a few days prior to the seismic events may be linked to the forthcoming events and most of them are nearly concurrent with atmospheric anomalies happening during the same day. This analysis also highlights that even in low-latitude areas it is possible to discern pre-earthquake ionospheric disturbances possibly linked with the imminent seismic events.


Introduction
Humanity has placed a huge effort to understand the generation mechanism of earthquakes even since ancient times with the ultimate goal to successfully predict them. This long-lasting endeavor which is extensively described by [1] begins with Anaximenes of Miletus and continues with Aristotle, who first conceived and described the "πνεṽµα" (pneuma, air in motion), an exhalation emanating from underground to the atmosphere and causing earthquakes. Other ancient Greek and Chinese philosophers also observed and reported pre-earthquake anomalies, such as fogs and clouds and the disappearance of water from wells. Similar pre-earthquake phenomena continued to be detected in the Renaissance and shortly later during the 16th, 17th and 18th century by a considerable number of scientists such as A. J. Buoni, Longobardo, Rossi, Mallet, Mercalli, however, it was only later at the beginning of 19th century when such anomalies were identified in order to verify the causal link between anomalous signals and earthquakes. By using TEC time series, [23] shown that the signal-to-noise ratio of a hypothetical TEC precursor signature is too low to be detected by their statistical analysis techniques, while the TEC analysis of [24] did not found any statistically significant difference in the spatial-temporal distribution of pre-earthquake ionospheric anomalies in the daytime and at night. However, the results of [23,24] as stated by them, do not disprove the potential existence of TEC earthquake precursors, but suggest the need for enhanced TEC ionospheric monitoring. Though [24,25] did not provide any clear evidence of statistically significant pre-earthquake TEC anomalies by using TEC time series deriving from Global Ionospheric Maps (GIM), they both highlight the fact that any possible preseismic or coseismic anomaly cannot be resolved by GIMs since their temporal and spatial resolution is relatively low and the final maps are interpolated both in time and space. Thus, short scale pre-earthquake anomalies cannot be detected.
Though the structure of ionospheric earthquake precursors is not well determined [23,24], in recent years a great progress on establishing the morphological features of these precursors has been made. Their structure is described by the LAIMC model. According to this model, all physical precursors of the LAIMC system are related to a chain of physical earthquake preparation processes as demonstrated by [2,20] which take place inside the earthquake preparation area defined as within cycle with a radius equal to R = 10 0.43M in km, M being the magnitude of the earthquake [26]. The underlying assumption of the LAIMC model is the increased radon emanation from underground to the earth's surface mainly up to 12 days prior to a seismic event. This radon release is capable of producing air ionization and water condensation nuclei by the produced ions, while the latent heat released from the condensation can modify atmospheric and thermal properties such as temperature, humidity and outgoing longwave infrared radiation (OLR). Changes of the air ionization and generation of large ion clusters affects the air conductivity of the planetary boundary layer resulting in the creation of an anomalous vertical electric field overactive tectonic faults, which through the global electric circuit between the earth's surface and the ionosphere, can cause abnormal variations of electron and ion plasma densities in the ionosphere. Such ionospheric irregularities are capable of penetrating to the magnetosphere through geomagnetic field lines due to their high conductivity.
There exists considerable experimental and theoretical evidence that supports the idea of increased radon exhalation prior to seismic events capable of inducing ionization of the air above the earth's surface. A list of about 180 radon earthquake precursory cases globally is reported in [27], whereas [28] examined pre-earthquake radon signals registered during the period 1983-2002 all over the globe too. Numerous studies of detected radon anomalies in ground water, soil and air since 1960s are also listed by [29,30]. Radon, a product of the uranium radioactive decay chain with a half-life of 3.8 days, is produced continuously and is always present in soil, rocks, surface water and any geological structure. During the earthquake preparation period the increased ground deformation causes high stress propagation in the soil that results in surges in radon emanation. Radon has the capacity to get transported for some distance, by diffusion, dissolution and transport by water, and convection, before it decays [27,29]. According to [31,32], a small variation in gas flow velocity into the ground causes important changes in radon emanation. The released radon is accumulated at the ground surface during nighttime owing to the low mixing height of the atmospheric boundary layer at night. Since the sudden radon accumulation results to a quick growth of ionization rate [33], it is anticipated that the nighttime radon storage will augment the ion conductivity [17,34].
The increased surface air ionization by the emanated radon can effectively modify the atmospheric electric field prior to a seismic event [35][36][37]. Experimental proof of this statement is the detection of electric field after an underground nuclear explosion [38] as well as the observed changes in electromagnetic properties within the earth-ionosphere guide after the emission of radioactive dust during the Chernobyl atomic plant accident [4]. Atmospheric electric field modifications associated with impending earthquakes have been reported by many authors [36,39,40].
In this context, the present work aims to investigate possible atmospheric and ionospheric anomalies preceding three large earthquakes (M > 6.6) that took place in Mexico during 2016 and 2017. A multi-technique analysis is performed including statistical and spectral analysis of the ionospheric Total Electron Content (TEC) and mapping of TEC and atmospheric chemical potential (ACP). To obtain increased statistical significance about our results we also examine all TEC anomalies during a long-term interval (the year 2017) enhancing also our present knowledge regarding the characteristics of seismo-ionospheric anomalies in particular over the low-latitude regions.

Materials and Methods
The relevant information about the earthquake characteristics was obtained from the earthquake catalog provided by the United States Geological Survey's (USGS) Earthquake Hazards Program and is presented in Table 1. The TEC values for the year 2017 as well as the period 1 to 22 January 2016 were retrieved from the International GNSS IGS network. We have processed RINEX files obtained from GNSS receivers according to the [41] algorithm, which computes vertical TEC (hereafter TEC) by using the relative phase delay between L1 and L2 GNSS frequencies while the signals broadcasted from each GNSS satellite propagate through the ionosphere. The output is given on a 30 s time resolution. In Figure 1 a map of the area under investigation is shown along with GNSS receiver stations, epicenters and the preparation zone of the three seismic events.
The methods we apply for the ionospheric and atmospheric analysis of possible preearthquake ionospheric anomalies are presented next. First, we try to identify possible significant TEC disturbances prior to the earthquake using a statistical method. Based on this method, if we assume a normal distribution with mean µ and standard deviation σ of TEC, we can create upper and lower bounds using µ ± 1.34σ. If the observed TEC exceeds the bounds for an epoch, then an abnormal signal is detected at a confidence level of 82% [42]. The method is applied for a period of 21 days before and one day after the event of 21 January 2016, whereas for the other two events in September 2017 we decided to apply the method for the whole year 2017 instead for a shorter period before each event and examine TEC anomalies prior to all the earthquakes with magnitude M ≥ 4.3 that occurred close to one GNSS station (GUAT) in Mexico within this year. The reason for this is to increase the reliability and obtain additional statistical significance of our results in regard to potential earthquake precursory phenomena as well as to acquire further information about their main features.
Ionospheric anomalies can also be caused by disturbed geomagnetic conditions. To exclude such anomalies, we carefully examine the geomagnetic bulletins provided by the National Geophysical Data Center (NGDC) (https://www.ngdc.noaa.gov/) as well as several geomagnetic and solar indices. More specifically, the Dst and AE geomagnetic indices provided by the World Data Center for Geomagnetism of Kyoto University (http://wdc.kugi.kyoto-u.ac.jp/wdc/Sec3.html) as well as the Solar index F10.7 and the component Bz of IMF obtained from the OMNI dataset of the Goddard Space Flight Center (https://omniweb.gsfc.nasa.gov/form/dx1.html ) are checked. Ionospheric anomalies can also be caused by disturbed geomagnetic conditions. T exclude such anomalies, we carefully examine the geomagnetic bulletins provided by th National Geophysical Data Center (NGDC) (https://www.ngdc.noaa.gov/) as well a several geomagnetic and solar indices. More specifically, the Dst and AE geomagneti indices provided by the World Data Center for Geomagnetism of Kyoto Universit (http://wdc.kugi.kyoto-u.ac.jp/wdc/Sec3.html) as well as the Solar index F10.7 and th component Bz of IMF obtained from the OMNI dataset of the Goddard Space Fligh Center (https://omniweb.gsfc.nasa.gov/form/dx1.html ) are checked.
Furthermore, the meteorological conditions for the periods under investigation ar considered by inspecting the meteorological synoptic maps so as to identify cyclones o weather fronts that may induce ionospheric disturbances [43] with similar characteristic to the possible earthquake associated ones. These maps are obtained from NOAA/ESRL Physical Sciences Division (PSD (https://psl.noaa.gov/cgi-bin/data/composites/printpage.pl).
Next, we apply the method of Ionospheric Precursor Mask (IPM) developed b Pulinets et al., (2002b) [44]. Considering that ionospheric variations show self-similarity, we construct a TEC pattern (IPM) for the study region to identify specifi spatial and temporal characteristics of TEC precursors for the certain region. For this w first calculate ΔTEC deviations from the current mean following the equation:  Furthermore, the meteorological conditions for the periods under investigation are considered by inspecting the meteorological synoptic maps so as to identify cyclones or weather fronts that may induce ionospheric disturbances [43] with similar characteristics to the possible earthquake associated ones. These maps are obtained from NOAA/ESRL, Physical Sciences Division (PSD) (https://psl.noaa.gov/cgi-bin/data/composites/printpage.pl).
Next, we apply the method of Ionospheric Precursor Mask (IPM) developed by Pulinets et al., (2002b) [44]. Considering that ionospheric variations show a self-similarity, we construct a TEC pattern (IPM) for the study region to identify specific spatial and temporal characteristics of TEC precursors for the certain region. For this we first calculate ∆TEC deviations from the current mean following the equation: where TEC is the current vertical TEC value at the fixed instant and TECa is the moving mean of vertical TEC assessed over 15 previous values of vertical TEC for the same instant. ∆TEC values are calculated for 15 days prior and 18 days after the largest earthquake in 2017 (8.2 Mw, 8th September) and are then shown in a plot (mask) where x-axis represents the number of days before and after the earthquake day (day 0) and y-axis corresponds to the local time. By repeating the process for a group of earthquakes occurring in the same region and then averaging ∆TEC datasets for each selected earthquake one may obtain a generalized ionospheric precursor pattern for typical earthquakes of the certain region. Detailed explanation of this method is given by [45].
Complementary to the statistical analysis of TEC time series and the TEC precursor pattern technique, the global ionospheric TEC maps (GIMs) are inspected so as to obtain a thorough picture about the spatial characteristics of pre-seismic ionospheric anomalies. GIMs are retrieved from the DRAWING-TEC project (Dense Regional And Worldwide International GNSS-TEC observation, http://seg-web.nict.go.jp/GPS/DRAWING-TEC/) and the International GNSS Service (ftp://cddis.gsfc.nasa.gov/pub/gps/products/ionex) where GIMs are made available in IONEX format.
The ionospheric analysis is completed by the application of Spectral analysis on slant TEC observations (following the satellite signal path) retrieved from several GNSS stations located inside the earthquake preparation zone of the three seismic events under investigation. The slant TEC was computed by using an algorithm developed by [46] which is based on equations provided by [47]. To exclude hardware biases from satellite and receiver hardware the method was performed on differential slant TEC data (dSTEC), taken as the difference of slant TEC measurements between two successive satellite epochs.
Spectral analysis is capable of detecting possible TEC disturbances occurring during geomagnetically disturbed periods, since by selecting the period T of TEC oscillations to be lower than 40 min, the geomagnetically induced TEC oscillations which are known to have periods of approximately 1 h are excluded. In addition to this, by using spectral analysis we are able to observe localized and short-scale TEC disturbances which appear inside the preparation area and even very close to the epicenter, as this analysis is performed on slant TEC measured simultaneously within several paths from each GNSS station to each satellite crossing the earthquake preparation zone. Apart from the location, the characteristics of the observed TEC oscillations are provided such as their period and amplitude of intensity.
The atmospheric analysis of earthquake precursory phenomena includes the investigation of Atmospheric Chemical Potential (ACP) in terms of time series on a certain station close to the epicenter and regional maps at particular pre-earthquake times. This is a new integral parameter introduced by [20] which comprises an indicator of radon variations prior to a seismic event. It is proved that latent heat for water molecules during phase transitions is equal to its chemical potential or to the work function when the molecule disconnects from water droplet [48]. The air ionization by surface emanated radon before an earthquake leads to the formation of additional ions in the atmosphere through nucleation and to latent heat release due to water condensation on ions and not only on water droplets. Therefore, prior to an earthquake, the total latent heat in the air needs to be corrected by adding an additional component due to the latent heat release from the newly created cluster ions from the radon ionization of the air. This is done by adding the ACP parameter (∆U) in the following equation of relative humidity H(t): where U(t) = Uo + ∆Ucos 2 t, with ∆U being the volume averaged correction of chemical potential, due to the external forcing of the environment, Uo is the chemical potential for pure water. Boyarchuk et al. 2006 [33] provided a mathematical relation for ∆U which is a function of air temperature at Earth's surface (T g ) and relative humidity of air (H): ∆U is the correction of total chemical potential of the water vapor in atmosphere, reflecting the formation of cluster ions. Increased ∆U values correspond to increased energy of water molecule attached with ions and to more stable and larger cluster ions. ∆U is greatly affected by weather conditions, topography and amounts of radon over a region. Nevertheless, temporal variability and spatial distribution of ∆U are undoubtedly linked to the earthquake preparation mechanism. In this work, ACP time series for a period of about 20 days before to 5 days after the three large seismic events of 2016 and 2017 in Mexico are analyzed. The values of ∆U can be derived from Equation (3) by employing NASA atmospheric assimilation models. It should be noted at this point that all the GNSS stations under investigation were located close to the epicenter, at a distance less than 700 km.

Results
In this section the results of the multi-technique and multi-parameter methodology to observe ionospheric and atmospheric perturbations possibly pertinent to the impending earthquake described earlier are presented and interpreted for the selected earthquakes that occurred recently in a low-latitude region.

Ionospheric Analysis
We examine abnormal ionospheric variations prior to the selected earthquakes in Mexico by applying statistical analysis on TEC as obtained from various GNSS stations located inside the earthquake preparation zone. For this, a long-term period of one year is considered. We then review global ionospheric TEC maps (GIMs) and we construct the Ionospheric Precursor Mask by utilizing TEC data from GIMs in order to understand salient temporal and spatial features of pre-earthquake ionospheric anomalies in particular over a low-latitude region where the ionosphere exhibits a peculiar behavior. We complete the analysis by performing spectral analysis on GNSS derived TEC observations which allow as to detect localized short-scale ionospheric anomalies prior to seismic events even during disturbed geomagnetic conditions.

Statistical Analysis on Ionospheric TEC
The diurnal TEC variations along with the corresponding upper, lower bounds (µ ± 1.34σ) and mean values are estimated for an interval of 33 days (from 20 August to 21 September 2017) for four stations (GUAT, SSIA, MANA, INEG) located inside the earthquake preparation zone and near the epicenter of both events that took place on 8 and 19 September 2017 ( Figure 2). Significant increase in TEC (close to or beyond the upper limit) is observed in all four GNSS stations on 3 September which lasted for about five hours and maximized at 21UT (16LT). Since the geomagnetic conditions were quiet during that day and since ACP anomaly and enhanced TEC wave oscillations close to the epicenter are also noticed during the same time (as it will be shown in later sections) we believe that the TEC anomaly on 3 September maybe linked with the forthcoming earthquake. This anomaly is clearly seen on the Ionospheric Precursor Mask ( Figure 3) which shows the percentage TEC differences from the running 15-day mean of the corresponding TEC values over the GNSS station UTON which is located close to the epicenter. In the same IPM one can identify large positive TEC deviations from the mean values also during 27 August and 7 to 8 September which are also seen in the TEC times series of all the rest GNSS stations (GUAT, SSIA, MANA, INEG), however they cannot be ascribed to the imminent 8.2 Mw earthquake on 8 September as unsettled (27 August) or disturbed (7-8 September) geomagnetic conditions were prevailing during these days. In particular, an intense geomagnetic storm was initiated on 7 September around eight hours before the earthquake and was followed by a second storm on 8 September.
Besides these large-scale TEC enhancements, a short-scale sharp increase of TEC is noted approximately three hours (around 2 UT) before the earthquake on 8 September in all four GNSS stations, which, though it occurred during a geomagnetic storm, can be considered to be linked with the impending earthquake, as geomagnetic storms do not prompt such short and sharp types of TEC anomalies as the one observed. To prove this postulation, in the next section we apply spectral analysis on TEC datasets by selecting the period of TEC wave oscillations to be less than 1 h so as to exclude geomagnetically induced TEC oscillations which are known to have periods larger than an hour. [49] also reported a short-lasting anomaly in magnetic field detected interestingly only by one of SWARM satellites (Charlie) 20 min before the main shock and exactly at the latitude and longitude of the same earthquake on 8 September. induced TEC oscillations which are known to have periods larger than an hour. [49] also reported a short-lasting anomaly in magnetic field detected interestingly only by one of SWARM satellites (Charlie) 20 min before the main shock and exactly at the latitude and longitude of the same earthquake on 8 September.  A second large earthquake (7.1 Mw) followed the 8.2 Mw earthquake only 11 days after (19 September) in the same region. Since the first 8.2 Mw event coincided with the occurrence of two successive intense geomagnetic storms, the second 7.1 Mw event 11 days later occurred within the long-lasting recovery phase of the storm and thus during unsettled geomagnetic conditions. Hence, the remarkable positive TEC enhancements which are found few days (14-15 September) and one day before this seismic event in all GNSS stations are most probably associated with the geomagnetic storms and not with the impending earthquake (Figures 2 and 4). On 18 September, one day prior to the earthquake, a large increase of TEC was observed in three GNSS stations located southeast of the epicenter. Though geomagnetic conditions were slightly unsettled during that day, the coexistence of ACP enhancements and of intensified TEC wave oscillations that will be presented in the next sections may point to a seismic origin of this anomaly. The inclination of the observed TEC anomalies southeast from the epicenter vertical projection complies with [4] who have found that in middle latitudes the ionospheric precursory phenomena are detected to the south (north) from the epicenter projection in the Northern (Southern) Hemisphere. A second large earthquake (7.1 Mw) followed the 8.2 Mw earthquake only 11 days after (19 September) in the same region. Since the first 8.2 Mw event coincided with the occurrence of two successive intense geomagnetic storms, the second 7.1 Mw event 11 days later occurred within the long-lasting recovery phase of the storm and thus during unsettled geomagnetic conditions. Hence, the remarkable positive TEC enhancements which are found few days (14-15 September) and one day before this seismic event in all GNSS stations are most probably associated with the geomagnetic storms and not with the impending earthquake (Figures 2 and 4). On 18 September, one day prior to the earthquake, a large increase of TEC was observed in three GNSS stations located southeast of the epicenter. Though geomagnetic conditions were slightly unsettled during that day, the coexistence of ACP enhancements and of intensified TEC wave oscillations that will be presented in the next sections may point to a seismic origin of this anomaly. The inclination of the observed TEC anomalies southeast from the epicenter vertical projection complies with [4] who have found that in middle latitudes the ionospheric precursory phenomena are detected to the south (north) from the epicenter projection in the Northern (Southern) Hemisphere.
The diurnal variations of TEC with its corresponding upper, lower bounds and mean values at stations INEG and GUAT, are evaluated for 18 days before, during and 2 days after the 6.6 Mw Mexico earthquake on 21 January 2016 (shown in Figure 4). A TEC increase beyond the upper bounds was detected 10 days before the earthquake, on 11 January at both stations and during approximately the same time period. The enhancement lasted for more than three hours at both stations (INEG inside and GUAT outside the earthquake preparation area) with a peak around 20 UT (15 LT), however, we cannot safely conclude the enhancements are due to the earthquake because unsettled geomagnetic conditions were prevailing at that time (Dst values were around −20 nT). In addition to this, such TEC enhancements can be related to the peculiarity of the ionosphere over Mexico which according to [50] is characterized by large positive TEC enhancements. Large positive TEC anomalies are also observed one day before the main shock on 20 January, which are attributed to the presence of a moderate geomagnetic storm.  The above statistical analysis of the ionospheric variations prior to the three earthquakes in Mexico clearly demonstrates the necessity of applying additional methods and examining more factors so as to derive concrete conclusions in regard to the detection of possible earthquake precursory phenomena and their characteristics. The fact that all three events coincided with the occurrence of geomagnetic storms and that the ionosphere over Mexico presents a distinct behavior, prompted us, for more reliance, to investigate TEC variations for a longer interval, the entire year 2017, so as to infer addi- The diurnal variations of TEC with its corresponding upper, lower bounds and mean values at stations INEG and GUAT, are evaluated for 18 days before, during and 2 days after the 6.6 Mw Mexico earthquake on 21 January 2016 (shown in Figure 4). A TEC increase beyond the upper bounds was detected 10 days before the earthquake, on 11 January at both stations and during approximately the same time period. The enhancement lasted for more than three hours at both stations (INEG inside and GUAT outside the earthquake preparation area) with a peak around 20 UT (15 LT), however, we cannot safely conclude the enhancements are due to the earthquake because unsettled geomagnetic conditions were prevailing at that time (Dst values were around −20 nT). In addition to this, such TEC enhancements can be related to the peculiarity of the ionosphere over Mexico which according to [50] is characterized by large positive TEC enhancements. Large positive TEC anomalies are also observed one day before the main shock on 20 January, which are attributed to the presence of a moderate geomagnetic storm.
The above statistical analysis of the ionospheric variations prior to the three earthquakes in Mexico clearly demonstrates the necessity of applying additional methods and examining more factors so as to derive concrete conclusions in regard to the detection of possible earthquake precursory phenomena and their characteristics. The fact that all three events coincided with the occurrence of geomagnetic storms and that the ionosphere over Mexico presents a distinct behavior, prompted us, for more reliance, to investigate TEC variations for a longer interval, the entire year 2017, so as to infer additional characteristics of seismo-ionospheric anomalies. In the analysis we include all the earthquakes with magnitude higher than 4.3 where the GUAT station is located inside their preparation zone. In general, anomalies are expected for earthquakes with magnitude greater than 5 [1,2], however, there is quite a number of studies which have reported anomalies prior to weaker earthquakes, such as [36] who identified a significant decrease of the atmospheric electric field during a 4.1 Mw earthquake in Portugal. Furthermore, in the one-year analysis in addition to ionospheric variations we examined the geomagnetic, solar and meteorological conditions in order to exclude TEC anomalies induced by such conditions. TEC variations along with the monthly mean and upper and lower confidence bounds are shown on a monthly basis in Figure 5.
Omitting the earthquakes occurring in the time frame of 11-21 May and April 2017 due to the missing TEC data, a total of 37 earthquakes (with M ≥ 4.3) which occurred in 2017 are selected for our analysis from the USGS catalog. First, to identify the days with disturbed geomagnetic conditions that may also be related to TEC anomalies, the Dst geomagnetic index timeseries is examined. Time frames having Dst value less than −10 nT are considered to be geomagnetically disturbed periods, which included not only magnetic storms but also unsettled geomagnetic conditions (having Dst between −10 to −20 nT). [51] shown that even a drop of Dst index to −15 nT is capable of inducing abnormal TEC variations. It is found that among the 37 examined earthquakes, 20 earthquakes exhibited pre-earthquake TEC anomalies during disturbed geomagnetic conditions, while in 6 earthquakes their respective observed TEC anomalies occurred during disturbed geomagnetic storms and close to storm conditions (grey boxes in Figure 5). Therefore, these anomalies were not regarded as precursors of impending earthquakes.
In addition to this, following [51][52][53][54], who postulate that the local TEC anomalies often reflect global changes of the ionization caused by the solar-driven 27-day variations during the rising or falling phase of solar activity, we also check the Solar index F10.7 (https://omniweb.gsfc.nasa.gov/form/dx1.html) during the examined year 2017 which is part of the falling phase of solar cycle 24. Furthermore, [51] also demonstrated that a drop of the IMF Bz down to −5 nT is sufficient to engender~20-30% excess of regional TEC over the 15-day median values. Therefore, we examine IMF Bz values during 2017 as well. AE index plots (https://omniweb.gsfc.nasa.gov/form/dx1.html) are also carefully inspected in order to identify any substorm which may also affect electron content in the ionosphere. Finally, taking into consideration the results in [48], which show the possibility of appearance of TEC anomalies prior to low-latitude earthquakes due to their closeness to the geomagnetic storms, we also isolated TEC anomalies which occurred close to geomagnetic storms.  The grey-shaded areas mark the days with geomagnetic unsettled or disturbed conditions, light blue areas mark TEC anomalies of seismic origin, green-shaded areas denote anomalies close to geomagnetic disturbed conditions, light-pink areas denote typical TEC enhancements over Mexico and light-yellow areas include TEC anomalies that could be either due to earthquake or to typical TEC enhancements over Mexico.
In this investigation, nine earthquakes presented pre-earthquake TEC anomalies around geomagnetic storms, as seen in Figure 5 (green box), out of which 4 earthquakes exhibited TEC anomalies not only around but also during these geomagnetic storms ( Table 2, column C). It is obvious that more than 65% of the examined seismic events take place during geomagnetic storms or around them. The comparison between the storm sudden commencement (SSC) time and earthquake occurrence time showed that most earthquakes happen during the recovery phase of a storm or during unsettled geomagnetic conditions. The 8.2 Mw earthquake occurred 8 h after the SSC time of the storm, whereas the 5 Mw earthquake on 31 August occurred exactly at the SSC time. Among the earthquakes we analyzed, few earthquakes occurred one day before or after the storm as well ( Table 2, column C).
In recent years, research has been initiated on the potential relationship between geomagnetic storms and earthquakes in the sense that solar activity can trigger earthquakes. During the 2020 EGU meeting, Pulinets and Khachikyan [19] pointed out that both cosmic rays and high energy particles precipitating from the radiation belt in times of strong geomagnetic storms may trigger earthquakes and showed some experimental results which supported this speculation. Additionally, [55] report that, undoubtedly, solar flares may act as a trigger for strong earthquakes. More studies on the propounded hypothesis that associates solar activity to earthquake activity needs to be carried out, and in the future, such research may add another dimension to the identification of pre-seismic ionospheric anomalies.
After the examination and exclusion of all the anomalies due to geomagnetic, solar and meteorological activity, it is found that four earthquakes showed only negative TEC anomalies (having TEC values less than the lower limit or much lower than the mean µ) that could probably be linked to the imminent earthquakes ( Table 2, column A). The absence of other physical phenomena (capable to cause such a TEC deviation) like tropical cyclones, volcano eruptions (occurred in 1-2 July), equatorial plasma bubbles (EPBs), solar eclipse, typical TEC enhancements in 20 • N over Mexico [50] during the pre-earthquake period also supported the possible linkage of these anomalies with these earthquakes. Recently, [7] confirmed the perturbing effect of volcanic and dust events on atmospheric and ionospheric properties. Further, in case of seismically induced ionospheric TEC anomalies, the anomaly duration is usually 4-6 h, but for very strong earthquakes the anomalies can last even up to 12 h as observed prior to the Good Friday earthquake of 9.2 Mw [1]. Here, significant negative anomalies lasting for few hours were observed on 12 and 13 June, and two and one days prior to the 6.9 Mw earthquake on 14 June. Past study of 1492 earthquakes during 1998 to 2014 with a magnitude greater than 5 has shown that sharper seismo-ionospheric disturbances are observed for earthquakes of Mw ≥ 6 [56]. Negative anomalies observed here prior to 4.6 Mw earthquakes might be of seismic origin as the epicenter of those earthquakes occurring on 15 January and 21 March are at a distance of only 41.56 km and 82.47 km respectively from the GUAT station. Similar results have been also shown by [57] in which TEC enhancements were observed prior to the 4.9 Mw Pipalkoti earthquake, as its epicenter was located very close to the stations ( Figure 5, light blue colored areas & Table 2, column A). It should be highlighted here that noone of the earthquakes in which possible earthquake-related ionospheric anomalies were detected, comprise aftershocks of a large preceding seismic event, as they occurred more than three days after large earthquakes. Therefore, the detected pre-earthquake ionospheric anomalies are most probably linked with these events and are not false precursors. To investigate in depth the negative TEC anomalies probably related to the above mentioned 4 earthquakes and elucidate their main features, we examine the global ionospheric TEC maps (http://stdb2.isee.nagoya-u.ac.jp/GPS/GPS-TEC/) having a high timeresolution of 10 min for all the days of the three months of 2017, January, March, and June, in conjunction with ionex maps (ftp://cddis.gsfc.nasa.gov/pub/gps/products/ionex). As the epicenters of the examined earthquakes are located at low-latitude (approximately ±20 • from the equator), the ionospheric equatorial anomaly is also expected prior to the main shock. Inspection of GIMs showed modifications in Equatorial Ionization Anomaly (EIA) structure several days before these four earthquakes which can be a harbinger of the forthcoming seismic events. Both crests of EIA disappeared and only one crest was formed over the equator on 13 January, the day when TEC negative anomaly was also detected. Such a "vanishing" of crest of EIA was also observed by [54] 10 days before the 7.8 Mw earthquake that occurred on April 16, 2016 in Ecuador. Several past studies have also reported similar equatorial ionospheric anomalies prior to low-latitude earthquakes [56][57][58][59][60][61]. Though the 4.6 Mw earthquake on 15 January is small to produce such an EIA identification, it should be noted that during the same time period larger earthquakes in the region took place such as the 5.8 Mw earthquake on 17 January. Figure 6 shows the EIA modification during afternoon local times (19)(20)(21)(22)(23)(24). Similar changes in the EIA structure were identified prior to the 4.6 Mw earthquake on 21 March. A number of large earthquakes of magnitude 5 to 6 also occurred during the same period (24-27 March), which could have induced the disappearance of the EIA crests, and the replacement by one crest closer to the equator on 15 and 18-20 March during the afternoon. It has been mentioned by [57] that such a disappearance of EIA crests is fundamentally a consequence of decrease in the vertical TEC. Finally, prior to the large earthquakes of magnitude 5 and 6.9 on 13 and 14 June, respectively, we found that the north crest of the EIA was enhanced largely on 1, 2, 8, 9, 15 June (geomagnetic quiet days) which corresponds to the TEC enhancements over the GUAT station. It was shown by [4] that changes in the global-scale ion density distribution are possible up to seven days before an earthquake event. As all the days mentioned above where EIA modifications were detected are geomagnetically quiet days, the observed modifications are not related to the storms. The negative TEC anomalies 1-2 days before the events (12-13 June) are most probably related to the imminent earthquake and not to solar flux variations, since they were observed only prior to earthquakes and not in any other geomagnetically quiet day of June.
Past studies have shown that TEC anomalies due to earthquake could be positive or negative [20]. The TEC time series of 2017 over Mexico region also shows prominent TEC intensifications on 5 and 10 June, 5 July and 31 August prior to seismic events. However, in this study, the peculiar characteristics of daily, seasonal and annual TEC variation over the Mexico region must be considered before deeming any positive anomalies as of seismic origin. In particular, [50] have concluded that short-lived TEC enhancement is one of the specific features observed in Mexico region, while [62] reported the existence of persistent positive TEC anomalies over the Gulf of Mexico. Therefore, even though the decreasing trend of solar flux in June ensured the enhancements on June were not due to solar activity, the observed positive anomaly still cannot be considered to be of seismic origin (Table 2, column B). Similarly, the TEC increase on 5 July and 31 August can be either due to the typical persisting TEC enhancements over Mexico or due to the temporal proximity with disturbed geomagnetic and solar flux conditions (Table 2, column D). As the physical mechanism which is responsible for such peculiar enhancements in Mexico region has not been identified yet, we are unable to discern the origin of such enhancements. Likewise, prior to 3 earthquakes that occurred on 28 February, 12 October and 4 November, large TEC enhancements occurred, sometimes for continuous days, lasting for about three or more hours, and cannot be concluded to be of seismic origin ( quake and not to solar flux variations, since they were observed only prior to earthquakes and not in any other geomagnetically quiet day of June. Past studies have shown that TEC anomalies due to earthquake could be positive or negative [20]. The TEC time series of 2017 over Mexico region also shows prominent TEC intensifications on 5 and 10 June, 5 July and 31 August prior to seismic events. However, in this study, the peculiar characteristics of daily, seasonal and annual TEC variation over the Mexico region must be considered before deeming any positive anomalies as of seismic origin. In particular, [50] have concluded that short-lived TEC enhancement is The one-year analysis of pre-seismic ionospheric TEC perturbations demonstrated that although most of the observed TEC anomalies are induced by many different sources such as geomagnetically disturbed conditions, solar variation, meteorological and other physical phenomena like for instance volcano eruptions and peculiar ionospheric TEC enhancements over the Mexico region, there are still TEC anomalies that cannot be linked to any of the above sources, but are most probably correlated with impending earthquakes. These present specific characteristics such as the modification of EIA over low-latitude regions, that may be either positive or negative and their duration can be between few minutes to few hours. Interestingly, a potential link between geomagnetic storm and earthquake occurrence was revealed; however, a deeper and more extensive investigation is needed to elucidate this hypothesis. It is further indicated from this analysis that additional attention is needed to discern seismically induced ionospheric anomalies from anomalies of other origin over equatorial and low-latitude regions, due to the particular irregularities that characterize the ionosphere over these regions.

Spectral Analysis of Ionospheric TEĊ
It has become evident that it is difficult to infer whether or not the observed TEC anomalies are correlated with the impending earthquake using statistical analysis only, even if one-year TEC variations were analyzed, not only because of the long-lasting disturbed geomagnetic conditions but also due to the peculiarity of the ionosphere over the Mexico region, that could drive TEC anomalies. Therefore, spectral analysis is performed on GNSS ionospheric measurements of TEC for a rationalized interpretation and confirmation of the potential earthquake-related TEC anomalies identified by the statistical analysis. The examined events are the three large earthquakes on 8 and 19 September 2017 (8.2 Mw and 7.1 Mw respectively) and 21 January 2016 (6.6 Mw). More specifically, spectral analysis is applied on slant TEC measurements obtained from GNSS stations located inside the earthquake preparation zones (GUAT, MANA, INEG) during the day of the earthquake and several days before (up to 12) as well as on seismically quiet days.
To further corroborate the observed TEC anomalies during 3 September (see Section 1) as of seismic origin, we carry out spectral analysis for TEC measurements on this day. The analysis clearly shows intensified short-scale TEC wave-like oscillations on 3 September and particularly at approximately the same time (around 21 UT) that large-scale TEC positive anomalies maximized. At that time also, the satellite path (PRN 25) is crossing over the region close and south of the epicenter at around 21 UT time (Figure 7). On the same day, an increase of ACP is observed which continues at 21 UT over the same area as shown in the next section. The approximate concurrence of both ionospheric and atmospheric perturbations prior to the main shock is supportive of the LAIMC model which describes that it is possible not only to detect a set of different earthquake precursory anomalies (atmospheric or ionospheric), but also precursors to appear in a specific time sequence with some delay as the height at which the recorded anomaly appears increases (first appear modifications in lower atmosphere, which will shortly later induce ionospheric changes in the upper atmosphere). This "arrow in time" is explained by [20] by presenting the example of the L' Aquilla earthquake on April 2009, where increased radon concentration was detected which was then followed by a change in surface air temperature, an OLR anomaly, and finally an ionospheric TEC disturbance. The sequence of anomalies prior to the seismic shock indicates a possible pre-earthquake process directed to a critical point that is the main shock and can be detected following a multi-parameter method. The location of the detected short-scale TEC anomaly south of the epicenter as mentioned earlier is characteristic of the low-latitude or equatorial earthquakes and is explained by [59,63].
To conclude whether the sharp, short-scale TEC enhancement that was noted a few hours before the 8.2 Mw earthquake on 8 September in all four GNSS stations (Figure 2) is linked to the event, a spectral analysis was performed. Intensified TEC wave fluctuations appear close to the area above the epicenter four hours before the event at around 1 UT (20LT) and continued up to approximately four hours after the earthquake having periods around T = 25 min (Figures 8 and 9). Similar TEC waves are observed in the spectral output of all three stations on 8 September and are also found on 7 September, one day before the earthquake, at the same local time (Figure 10). This self-similarity and localtime dependence of seismo-ionospheric variations is reported by [64] and has facilitated the development of a new pattern recognition method for the detection of ionospheric anomalies, named precursor mask for earthquakes as explained at the previous section. The TEC wave fluctuations obtained from spectral analysis on 7 and 8 September are also compared with the corresponding spectral output during "quiet" days 27 August and 6 September 2017 when seismic activity was absent or significantly reduced. Unlike, 7 and 8 September no such TEC fluctuations are found after 2 UT (21 LT) during these days ( Figure 11) and thus it can be concluded that they are most probably related to the 8.2 Mw earthquake on 8 September. This is also confirmed by the investigation of the pre-earthquake atmospheric anomalies shown in the next section.
Geosciences 2021, 11, x FOR PEER REVIEW 17 of 29 ionospheric changes in the upper atmosphere). This "arrow in time" is explained by [20] by presenting the example of the L' Aquilla earthquake on April 2009, where increased radon concentration was detected which was then followed by a change in surface air temperature, an OLR anomaly, and finally an ionospheric TEC disturbance. The sequence of anomalies prior to the seismic shock indicates a possible pre-earthquake process directed to a critical point that is the main shock and can be detected following a multi-parameter method. The location of the detected short-scale TEC anomaly south of the epicenter as mentioned earlier is characteristic of the low-latitude or equatorial earthquakes and is explained by [59,63]. To conclude whether the sharp, short-scale TEC enhancement that was noted a few hours before the 8.2 Mw earthquake on 8 September in all four GNSS stations ( Figure 2) is linked to the event, a spectral analysis was performed. Intensified TEC wave fluctuations appear close to the area above the epicenter four hours before the event at around 1 UT (20LT) and continued up to approximately four hours after the earthquake having periods around T = 25 min (Figures 8 and 9). Similar TEC waves are observed in the spectral output of all three stations on 8 September and are also found on 7 September, one day before the earthquake, at the same local time ( Figure 10). This self-similarity and local-time dependence of seismo-ionospheric variations is reported by [64] and has facilitated the development of a new pattern recognition method for the detection of ionospheric anomalies, named precursor mask for earthquakes as explained at the previous section. The TEC wave fluctuations obtained from spectral analysis on 7 and 8 September are also compared with the corresponding spectral output during "quiet" days 27 August and 6 September 2017 when seismic activity was absent or significantly reduced. Unlike, 7 and 8 September no such TEC fluctuations are found after 2 UT (21 LT) during these days ( Figure 11) and thus it can be concluded that they are most probably related to the 8.2 Mw earthquake on 8 September. This is also confirmed by the investigation of the pre-earthquake atmospheric anomalies shown in the next section.         Like the 8.2 Mw earthquake on 8 September, the 7.1 Mw event on 19 September, which was the second largest earthquake in 2017, also occurred during disturbed, though less strong, geomagnetic conditions and thus, it was impossible to detect any large-scale TEC anomaly by employing only statistical analysis. The application of spectral analysis, though on slant TEC observations deriving from GNSS stations located inside the earthquake preparation zone reveals the existence of intensified TEC wave fluctuations over the greater epicenter area during five and one days, as well as one hour before the event at around 17 UT (12 LT) ( Figure 12). These fluctuations occurred at the same time every day, between 16-17 UT (11-12 LT) and around 18-23 UT (13)(14)(15)(16)(17)(18) and their periods are around 25 min. These periods are also corresponding to the periods of IGWs, which according to the numerical simulations of [21] can induce a positive effect on TEC over the epicenter area. Since the observed TEC anomalies in our analysis are stationary and are lasting for few hours, we cannot derive unambiguous conclusions concerning the IGWs effect on the ionosphere. Additional investigation would be greatly desirable in order to come up with more solid conclusions.
The spectral results in conjunction with the atmospheric anomalies analysis presented later allow us to consider that the observed anomalies are most probably of seismic origin.
The large 6.6 Mw earthquake on 21 January 2016 also took place during a large geomagnetic storm and therefore, a multi-method approach was also adopted so as to detect possible ionospheric precursors. The inspection of spectrograms during the earthquake day indicated enhanced TEC wave fluctuations above the area close to the epicenter one hour prior to the event that is on 17 UT (12 LT), which continue up to three hours after the earthquake. In the same area, intensified TEC wave fluctuations appear during the interval 6-9 UT (1-4 LT) in several days before the earthquake (11,17,19 and 20 January), indicating that potential seismo-ionospheric precursors occur in certain time intervals during the day. According to [4], this is characteristic for a given place and for different seismic zones ( Figure 13). Considering finally the observed atmospheric anomalies prior to the event that are described analytically in the following section, it can be stated that these TEC disturbances may be linked to the 6.6 Mw earthquake on 21 January 2016.   In this point we should mention that for all three examined earthquakes intensified and sharp TEC oscillations are detected, which occur on the same time interval every day around sunset and sunrise times. They have recurrence periods similar to the seismically induced wave-like TEC fluctuations, however, they are sharp and short-scale variations unlike the ones of seismic origin which have spherical wave-like periodic structure. We attribute these fluctuations to the solar terminator transition which is known to be a source of TEC undulations with periods between five minutes to one hour and amplitudes ranging from 0.05-0.1 TECU [65,66]. Furthermore, TEC wave-like fluctuations which have similar characteristics with that of seismic origin can be induced by strong meteorological phenomena such as tropical cyclones and weather fronts [43]. For this reason, we have checked the weather maps for all the earthquake cases under investigation to exclude such TEC wave fluctuations.

Atmospheric Analysis
Atmospheric chemical potential (ACP) is a parameter which presents temporal and spatial variation within the time interval of 1 month to few days before strong earthquakes [45]. Along with statistical analysis of vertical TEC and spectral analysis on slant TEC, investigation of ACP time series on stations in the vicinity of epicenters and of ACP regional maps may support the identified anomalies in the ionosphere as possible earthquake precursors. ACP indicates how much the latent heat in the boundary layer of the atmosphere changes in the presence of an additional source which is the ionization of air by emanated radon prior to an earthquake. According to the LAIMC model, this radon release from active faults leads to ionization of air which contributes significantly to processes of energy transformation in the atmosphere [18]. The ions formed in the atmosphere undergo ion hydration (also called ion-induced nucleation, IIN) in which the molecules of water vapor and the ions join. The process results into the formation of large, charged cluster ions which become the centers of condensation [45]. The chemical potential of pure water molecules during condensation is equal to the latent heat, but the chemical potential of water vapor molecules attached to the ions is different than that of pure water drops. Taking into account this discrepancy in chemical potential, [20] introduced ACP which is a correction of the chemical potential of the water vapor in the atmosphere representing latent heat due to IIN. The release of latent heat due to vapor condensation on ions justifies the increase in air temperature and decrease in air humidity observed prior to earthquakes [20,67]. Increase in air temperature and drop in relative humidity might be detected over the earthquake preparation area [48] before strong earthquakes in Mexico, as the Colima earthquake 2003, the Oaxaca earthquake 1999 and the Michoacán earthquake 1985 [68].
ACP measurements in Kamchatka peninsula were analyzed by [69], considering 64 Kamchatka earthquakes > 6 Mw occurring in a period of more than 10 years. It was observed that the ACP maximum occurred within a month before the main shock and the minimum ACP was seen few days before or few days after the day of earthquake. Since the temporal dynamics of ACP observed prior and after the earthquake resemble radon variation, [69] proposed ACP measurements as a good substitute of radon measurement. Evidence of such temporal variation of ACP was observed prior to various earthquakes like the Van 2011 earthquake, the Virgina 2011 earthquake, the Northwestern Iran 2012 earthquake [70], and the Lushan 2013 earthquake [20].
In this context, and to further verify the observed TEC anomalies as of possible seismic origin, we also analyze here the ACP time series and maps prior to the three selected major earthquakes in Mexico. Figure 14a shows an ACP maximum over the point located at 16.02 • N, 266.13 • E and close to the epicenter on 3 September, which is five days before the main shock. On the same day increased ACP was observed over the epicentral region of the 8.2 Mw earthquake which continued at 21 UT of the same day (Figure 14b). The ACP maximum is synchronous with the large and short scale TEC anomalies observed at 21UT in the TEC time series and spectrogram. The same time series of ACP in Figure 14a also shows significant drop on 10 September 2017 (2 days after the 8.2 Mw earthquake). thermore, multi-parameter approach has been employed in [2] for analysis of the 6 Mw Napa Earthquake of 2014, the 6 Mw Taiwan and the 7 Mw Kumamto earthquakes o 2016, where along with ACP measurements, analysis of OLR and differential TEC (dTEC have also been carried out. Likewise, [72] have also used multi-method approach for th 7.9 Mw Wenchuan and the 7.3 Mw Yutian earthquakes in 2008. Through these and othe studies with similar methodological approach, it is demonstrated that pre-earthquak transient anomalies (OLR, ACP, dTEC) appear primarily with a short time-lag from hours up to a few days and lead to a critical point of a pre-earthquake mechanism whic is the seismic shock. The sequence of anomalous pre-earthquake phenomena was eviden also from our analysis, as mentioned earlier.  In Figure 14c two spikes are obvious, eight and seven days before the 7.1 Mw earthquake on 19 September (i.e., on 11 and 12 September respectively). This could be related to the impending earthquake. Similar ACP anomaly was observed six days before the Van 2011 earthquake [70]. The spatial distribution of ACP over the epicentral region at 3 UT, one day before the earthquake (18 September 2017) shows enhancements very close to the epicenter, especially towards the west (Figure 14d). Unfortunately, a lot TEC data for 11 and 12 September were missing and thus it was not possible to run spectral analysis and check if there is concurrency between ACP and TEC anomalies on those days.
In Figure 14e very high values of ACP on 16 and 17 January, 2016 are observed, which took place five and four days before the seismic event of 21 January and could be of seismic origin. In addition, the spatial variation of ACP at 9UT on 16 January showed a large ACP increase in the north-west region from the epicenter (Figure 14f). At the same time, the TEC spectrogram on 16 and 17 January demonstrated enhanced TEC fluctuations during 4-7 UT and 14-18 UT over the epicentral area.
In our analysis, ACP maxima were noted at different periods but within 4-8 days before the seismic event. The frequency distribution of the time of the main maximum ACP by [69] showed that most of the earthquakes among the 64 Kamchatka earthquakes exhibited an anomalous maximum four days before the main shock. Some past studies have reported ACP anomalies even earlier, as in [71], where steady increase in ACP was observed 12 days before the 6.2 M Central Italy earthquake on 24 August 2016. Furthermore, multiparameter approach has been employed in [2] for analysis of the 6 Mw Napa Earthquake of 2014, the 6 Mw Taiwan and the 7 Mw Kumamto earthquakes of 2016, where along with ACP measurements, analysis of OLR and differential TEC (dTEC) have also been carried out. Likewise, [72] have also used multi-method approach for the 7.9 Mw Wenchuan and the 7.3 Mw Yutian earthquakes in 2008. Through these and other studies with similar methodological approach, it is demonstrated that pre-earthquake transient anomalies (OLR, ACP, dTEC) appear primarily with a short time-lag from hours up to a few days and lead to a critical point of a pre-earthquake mechanism which is the seismic shock. The sequence of anomalous pre-earthquake phenomena was evident also from our analysis, as mentioned earlier.

Synopsis and Conclusions
An important concern to observe pre-earthquake phenomena has been declared even since ancient times in many different scientific fields independently, while in the last decades a vast endeavor has been invested to investigate such precursory phenomena by combining observations and knowledge of different physical characteristics describing several parts of the earth-atmosphere system (mainly lithosphere, troposphere, ionosphere and magnetosphere) considering that such phenomena comprise part of a unified earthquake preparation process initiated from underground and evolving within the atmosphere and up to magnetosphere several days before the main seismic shock. A huge number of signals preceding large earthquakes has been reported which are of diverse origin.
On this basis, we aimed at analyzing simultaneously the ionospheric (TEC) and atmospheric (ACP) abnormal variations prior to three large earthquakes in Mexico, by applying different methods (statistical and spectral analysis on TEC, GIMs, IPM, time series and regional maps of ACP) so as to provide additional evidence and gain new insights of the pre-earthquake physical mechanism that leads to earthquake precursory phenomena. In addition, TEC measurements during a long interval (one year) were examined which allowed us to compare all pre-earthquake ionospheric anomalies during 2017, identify their different sources and features and hence isolate the possible seismically-induced anomalies with high confidence. All factors which may cause ionospheric anomalies were taken into consideration for this, including solar and geomagnetic variations, meteorological and other physical phenomena such as solar terminator passages, volcano eruption, and typical TEC enhancements over Mexico.
The results of this analysis imply that it is feasible to capture potential earthquakerelated ionospheric and atmospheric anomalies by following a multi parameter and multimethod approach, despite the fact that all three selected earthquakes coincided with the occurrence of geomagnetic storms and that long-lasting disturbed geomagnetic conditions were prevailing within the entire year 2017. In this respect, both large and short scale ionospheric TEC perturbations are found to occur few hours (1-3) up to few days (10 days) prior to the examined earthquakes. With the aid of spectral analysis on TEC it is shown that short scale anomalies are located in the vicinity of the epicenter and shifted primarily southeast of its vertical projection. They also present periodic structure and periods around 20-25 min. In most cases TEC disturbances are accompanied with ACP enhancements over the same area. This fact may be considered as an experimental evidence of the evolution of the earthquake preparation stages of a unified system towards the main shock critical point. The observed pre-earthquake TEC anomalies during the whole year 2017 (of total 37 earthquakes) are predominantly positive, although prominent negative TEC anomalies were also detected. As discussed, such TEC depletions are most probably linked to the modification of the Equatorial Anomaly structure several days before the events, a typical characteristic of equatorial and low-latitude earthquakes. The ionosphere over these regions exhibits large irregularities as it is shown also in the present study, and therefore particular attention should be placed to discern earthquake related ionospheric anomalies from those induced by different physical phenomena over these regions.
Through the current investigation, the need to establish a network of different stations and instruments such as GNSS receivers, radon, temperature, relative humidity and atmospheric electric field sensors in a seismically active region, is exemplified. This scenario will facilitate a systematic and uninterruptable monitoring that will provide the data necessary to obtain a comprehensive picture of the underlying physical mechanism of the earthquake preparation period. Funding: Part of this work was funded by COST Action "Atmospheric Electricity Network: coupling with the Earth System, climate and biological systems (ElectoNet)", supported by COST (Cooperation in Science and Technology). Part of this work was also funded by Frederick University.
Institutional Review Board Statement: Not applicable.