Pre-Seismic Irregularities during the 2020 Samos (Greece) Earthquake (M = 6.9) as Investigated from Multi-Parameter Approach by Ground and Space-Based Techniques

Abstract

We present a comprehensive analysis of pre-seismic anomalies as computed from the ground and space-based techniques during the recent Samos earthquake in Greece on 30 October 2020, with a magnitude M = 6.9. We proceed with a multi-parametric approach where pre-seismic irregularities are investigated in the stratosphere, ionosphere, and magnetosphere. We use the convenient methods of acoustics and electromagnetic channels of the Lithosphere–Atmosphere–Ionosphere-Coupling (LAIC) mechanism by investigating the Atmospheric Gravity Wave (AGW), magnetic field, electron density, Total Electron Content (TEC), and the energetic particle precipitation in the inner radiation belt. We incorporate two ground-based IGS GPS stations DYNG (Greece) and IZMI (Turkey) for computing the TEC and observed a significant enhancement in daily TEC variation around one week before the earthquake. For the space-based observation, we use multiple parameters as recorded from Low Earth Orbit (LEO) satellites. For the AGW, we use the SABER/TIMED satellite data and compute the potential energy of stratospheric AGW by using the atmospheric temperature profile. It is found that the maximum potential energy of such AGW is observed around six days before the earthquake. Similar AGW is also observed by the method of wavelet analysis in the fluctuation in TEC values. We observe significant energetic particle precipitation in the inner radiation belt over the earthquake epicenter due to the conventional concept of an ionospheric-magnetospheric coupling mechanism by using an NOAA satellite. We first eliminate the particle count rate (CR) due to possible geomagnetic storms and South Atlantic Anomaly (SAA) by the proper choice of magnetic field B values. After the removal of the statistical background CRs, we observe a significant enhancement of CR four and ten days before the mainshock. We use Swarm satellite outcomes to check the magnetic field and electron density profile over a region of earthquake preparation. We observe a significant enhancement in electron density one day before the earthquake. The parameters studied here show an overall pre-seismic anomaly from a duration of ten days to one day before the earthquake.

1. Introduction

The seismic hazards and their preparation mechanism are extremely complex. The overall mechanism and its possible outcomes not only have huge outbursts of mechanical energy but have a wide range of physical and chemical processes attached to the lithosphere, atmosphere, and ionosphere [1]. Several studies have reported that the processes of preparation for earthquakes (EQs) are confined within a preparation or critical zone [2,3]. Electro-kinetic phenomena, emission of radioactive particles, thermal irregularities, emissions of electromagnetic signals, and many other physical processes have been identified as the short-term precursory phenomena of impending EQs. Numerous studies have already established the pre-, co-, and post-seismic anomalies using different parameters [1,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. It is found that the pre-EQ processes are nonlinear, anisotropic, and can have a dependency of multiple parameters. Based on such an idea, a coupling mechanism was established to understand the pre-seismic process that begins well before the mainshock known as the Lithosphere–Atmosphere–Ionosphere coupling (LAIC) mechanism [1]. This coupling (LAIC) mechanism is detected through various channels such as thermal [1], acoustic [21,22], and electromagnetic [23]. Many ground- and space-based instruments are used to detect these disturbances [1,24,25]. The anomalies are observed using ground-based observations with the Very-Low-Frequency (VLF)/Low frequency (LF) radio wave, Critical frequency of F2 layers ($f0F_2$), Total Electron Content (TEC), Ultra Low-Frequency (ULF) emissions, etc. [26,27,28,29,30,31,32,33,34]. In parallel, a lot of space-based techniques have been used to investigate the LAIC mechanisms. A lot of satellite observations provide the atmospheric ion-concentration, electron density, magnetic field components, etc. Using Detection of Electromagnetic Emissions Transmitted from EQ Regions (DEMETER) and CHAllenging Minisatellite Payload (CHAMP) satellite data, Ryu et al. [35] have computed electron temperature, electron concentration, ion concentration, and ion temperature to study the perturbation in the ionosphere due to the Honshu EQ which happened in Japan on September 2004 and have reported the ion temperature enhancement one week before the mainshock. Ryu et al. [36] showed the Equatorial Ionospheric Anomaly (EIA) one month before the Wenchuan EQ (12 May 2008) using DEMETER satellite data. They have concluded that this kind of anomaly happened due to the external electric field generated from the EQ epicenter that perturbs the E × B drift which already exists. There are several reports on an investigation of pre-seismic anomalies based on multi-parametric approaches in recent years. Pulinets et al. [37] have used the GPS (Global Positioning System)-TEC, Outgoing Longwave Infrared Radiation (OLR), Atmospheric Chemical Potential (ACP), and the b value from the Gutenberg–Richter Frequency Magnitude Relationship (FMR) simultaneously. They have observed depletion in the cross-correlation coefficient in GPS-TEC one day before, OLR anomaly 19 days before, anomalies in ACP distribution 10 days before, and decrement in b value for the Kamchatka EQ on 20 March 2016. In De Santis et al. [38], they have used four parameters, specifically skin temperature, methane, total column water vapor, and aerosol optical thickness. They observed positive increment in all of these chosen parameters before the California EQ on 5 July 2019. They have also observed the disturbance related to this EQ in the Swarm satellite-A on 3 June 2019, in the Y component of magnetic field variation. This anomaly was observed in the nighttime period. They have also observed the disturbances in the electron density profile from the same track. In Chetia et al. [39], they have reported anomalies in three parameters to study the ionospheric perturbation related to the Kokrajhar, Assam EQ on 12 September 2018. The examined parameters were TEC, geomagnetic total field intensity, and soil radon. They noticed significant anomalies in all parameters before the corresponding EQ. Recently, Piersanti et al. [40] have shown an extensive research on this multiparametric study using Magnetospheric-Ionospheric-Lithospheric Coupling MILC) model and justified their various parameters for 2018 Bayan EQ. In order to make an accurate comparison between the co-seismic and pre-seismic observations, they first separated the data set based on relative to the co-seismic and pre-seismic observations and then analyzed each of the four parameters: atmospheric oscillations (AGW), ionospheric plasma (TEC), electric field perturbations, and magnetospheric FLR eigenfrequencies. In this manuscript, we use ground-and space-based observation of some of the well-established parameters such as TEC, Atmospheric Gravity Waves (AGW), energetic particle bursts in the inner radiation belt, magnetic field intensity, and ionospheric electron density profile during the 2020 Samos earthquake (EQ) (M = 6.9, epicenter (37.9001${}^{\circ }$ N, 26.8057${}^{\circ }$ E)in Greece on 30 October 2020.

The ionospheric disturbances can be observed by using a dual-frequency GPS receiver to understand the behavior of TEC in the ionosphere. The total number of free thermal electrons present between the GPS satellite and the receiver is known as TEC. It is also known as columnar electron density. As the ionosphere is a dispersive medium for electromagnetic waves, it introduces a time delay in radio signals during the propagation of GPS satellite signals [41], which can be determined by comparing the two GPS satellite signal frequencies. Absolute TEC values were calculated using dual-frequency GPS receivers in this study. Numerous researchers have already established that TEC can be used as one of the seismic precursor parameters. The electron density variation before 1999 M = 7.6 Chi-chi quake is first reported by Liu et al. [27] utilizing TEC processed from the GPS receiver. They found that there was a decrement of TEC during the evening time frame on days 1, 3, and 4 preceding this EQ. In this paper, we utilized the method described in [28]. They confirm this outcome through a statistical analysis of a global ionosphere map (GIM) TEC during the 20 M ≥ 6.0 EQs in Taiwan from September 1999 to December 2002. After that, more comparable work has been reported by utilizing GIM to examine TEC peculiarities before huge EQs, and numerous conspicuous outcomes have been found. At that point, the factual investigation before 17 M ≥ 6.3 seismic EQs during the long-term time frame from 1 May 1998 to 30 April 2008 appeared from Liu et al. [42] utilizing GIM TEC. They tracked down that the TEC is diminished 3–5 days before the EQ over the focal point. The seismo-ionospheric antecedents of odd abatements in TEC, which seem to be five days preceding 26 December 2004 M9.3 Sumatra–Andaman EQ are accounted for in [43]. There is additionally some amazing perception around both focal point and its form point during the 2008 M 8.0 Wenchuan EQ [44,45,46,47]. Akhoonzadeh et al. [48] revealed that the measurable outcomes show the positive and negative oddities in both of DEMETER and TEC information during days 1–5 preceding all determined EQs during calm geomagnetic conditions are exceptionally viewed as seismo-ionospheric forerunners. The TEC over the focal point or form point accomplish their most extreme qualities on the day before the EQ over the mid scope area, while the northern peak of tropical ionization peculiarity (EIA) moves poleward [49]. There are various prompt perceptions of seismo-ionospheric inconsistencies by utilizing GPS-TEC in situ electron density from circling satellites [50,51,52]. There is an enhancement of TEC over the epicenter 1 day before 12 January 2010, M = 7 Haiti EQ. Tao et al. [53] contemplated the ionospheric varieties over the focal point of 17 July 2006 M = 7.7 south of Java seismic EQ. The outcomes show that seismo-ionospheric inconsistencies in the GPS-TEC and in situ plasma density happen at practically similar times over the epicenter by utilizing GPS-TEC and DEMETER plasma information. Vita et al. [54] showed the identification of ionospheric GPS TEC anomalies before the EQ in Sumatra between 2007–2012 using the correlation technique. Moreover, there are some additional works on TEC variation during the EQ in the Indian-subcontinent region. Kumar et al. [55] examined the ionospheric TEC because of the seismic EQ of Tamenglong on 3 January 2016. They contemplated the Tamenglong EQ (M = 6.7, 3 January 2016) by utilizing the information from the stations at Lhasa, China, Hyderabad, India (17.410${}^{\circ }$ N, 78.55${}^{\circ }$ E), and Patumwan, Thailand (13.730${}^{\circ }$ N, 100.53${}^{\circ }$ E). They tracked down huge augmentation during the five days before the EQ and decrement of TEC after the mainshock of the EQ. Sharma et al. [56] proposed an ionospheric TEC demonstrating to notice the variety of TEC during the time of 0–8 days before EQs: 1 April 2015, Pipalkoti EQ (M = 4.9, 25 April 2015) Nepal seismic EQ (M = 7.8), April 26, 2015 (M = 6.7) and 12 May 2015 (M = 7.3) Nepal EQ. In addition, they also tracked down that the estimation of TEC was lower before 13–14 days from the initial two seismic EQs. Generally speaking, the examination has uncovered that low TEC that followed a few high TEC are well associated with the seismic events in the Himalayan region.

The traveling ionospheric disturbances (TIDs) are created by EQs and tidal waves. Komjathy et al. [57] had identified “normal” TIDs across Japan around 5 h following the 2011 Tohoku EQ. The variety of observed TEC irritations can be predicted using JPL’s ongoing Global Assimilative Ionospheric Model (GAIM) system, which includes wave-created gravity waves, auroral rotation, regular TIDs, and tropical fluctuations. The variety of observed TEC irritations can be predicted using JPL’s ongoing Global Assimilative Ionospheric Model (GAIM) system, which includes wave-created gravitational waves, auroral rotation, regular TIDs, and tropical fluctuations. Amin et al. [58] have reported that lightning-induced AGW can be easily detected from the GPS-TEC by utilizing a filtration method on it. After computation of the fluctuations, the spectral analysis has given the wave-like structures. He studied the lightning of South Africa during 2012 and noticed the connection between ionospheric abnormalities and surprising lightning occasions, and this phenomenon might generate the gravity wave that travels up to the F region and modulated ionospheric electron density. Oikonamou et al. [59] made a broad examination of ionospheric TEC antecedents identified with the M = 7.8 Nepal and M = 8.3 Chile EQs in 2015 based on spectral and statistical analysis and discover the wave-like structures before the EQ. It was shown that peculiar TEC patterns appeared from a few days to 2–3 h prior to the events that lasted up to 8 h, but an extended TEC wave-like pattern with periods of 20 or 2–5 min was also recognized, which could be correlated with the approaching seismic events.

The formation of the AGW is furthermore expected as a fundamental part of phenomena antecedent of seismic events. The fundamental acting agent of the acoustic channel of LAIC is AGW, which can show up on account of the atmospheric oscillation close to the epicentral zone of the relating shudder, this effect going to the upward heading and irritating the ionosphere. Different front-facing frameworks and wind currents are the primary justification for making Gravity Waves (GWs). These GWs are the fundamental system to move the energy from the lower atmosphere to the upper stratosphere and mesosphere. An AGW is considered as quite possibly the main wave in the climate on account of its solid effect on nearby and worldwide environmental constructions. These motions could be produced from the difference in ground movements, temperature, and pressing factors. Ground-based estimations have been completed to acquire wind and temperature information to look at AGW activity [60]. AGW is produced a few days before the EQ and engenders along the upward bearing [14,61,62]. During EQ preparation, varieties in temperature, conductivity, and pressing factor bring about AGW generation in the atmosphere.

The speculation of AGW excitation due to seismic events before an EQ has been reported by [63,64,65]. Murayama et al. [66] has observed that the variation of AGW energies are directly related to the jet stream by using the middle and upper environment (MU) radar observational data in Japan. AGW development and its relationship with convections during Indian southwest tempest were studied by [67,68] using the Mesosphere–Stratosphere–Troposphere (MST) radar at Gadanki, India. Tsuda et al. [69] used a GPS-temperature profile to investigate the overall allocation of potential energy over mid-scope and showed that potential energy is more imperative in winter seasons. Korepanov et al. [14], with the help of surface ecological squeezing factor and appealing field data, assumed that AGW can be a huge limit in the seismo-ionospheric study. Zhang et al. [70] and Yang et al. [71] used a Sounding of the Atmosphere utilizing Broadband Emission Radiometry (SABER) instrument installed on the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) satellite temperature profile to examine $E_p$ related to AGW. Nakamura et al. [72] have played out a similar assessment and endeavored to find the relating seismogenic sway for certain EQs. For the 2004 Niigata–Chuetsu EQ (M = 6.8), wavelet examination of those limits showed changes periods of 10 to 100 min which is in the extent of AGW. Yang et al. [21] have first announced that the AGW theory can be utilized as EQ antecedent by utilizing the ERA5 temperature profile for the 2016 Kumamoto EQ. They noticed the abnormalities in the AGW movement only 4–6 days before the EQ. After that, in Yang et al. [22], they have detailed similar speculation for an alternate oceanic 2011 Tohoku seismic utilizing both the SABER and ERA5 temperature profiles and contrasted the outcomes against the Kumamoto EQ ones. In our past examination [73], we have detailed that the geomagnetic storm that happened around the EQ can not contaminate the AGW processed from the SABER temperature profile and checked our outcome with the AGW acquired from the VLF signal.

Past studies reflect that the perturbation in the ionospheric and magnetospheric height is induced by seismic events. This inducing perturbation occurred within the EQ preparation zone before the seismic event [33]. VLF electromagnetic discharges are known to be created a couple of hours before the large EQs in their epicenter region or above them [74]. When it travels through the atmosphere and the ionosphere, it might disturb the high energetic particles traveling path in the Van Allen Radiation Belt (VAB). This makes the high energetic particles to come down to the height of the Low Earth Orbital (LEO) satellite which is below VAB. Now, if the mirror points of the particles are higher than 100 km, then the high energy particles drift across the magnetic field lines (L-shell) of the earth [75]. This unusual disturbance consists of the entire perturbed L-shell and the time frame of particle longitudinal drift. At the point when an LEO satellite crosses the perturbed L-shell, an onboard particles detector detects a sharp enhancement in particle count numbers. For this situation, the L-shell of particle bursts agrees with the L-parameter of the local disturbance in the radiation belt. Thus, the investigation of particle bursts gives the chance to determine the zone of the radiation belt particle disturbance and subsequently to look for the generation of a seismic event [76,77,78]. Several attempts are taken to investigate the high energetic particle number enhancement due to seismic events. In 1980, the first proof of EQ-induced PBs is given by [79,80]. They showed short-term enhancement in the particle counts number in the inner radiation belt due to EQs. The conclusion was made based on the result using MARIA-1 experiment [81].

Fidani et al. [82] and Battiston et al. [83] presented the results of particle fluctuation related to a seismic event. For this detection of particle count number, they have used National Oceanic and Atmospheric Administration (NOAA) satellite data. The NOAA and the National Aeronautics and Space Administration (NASA) mutually built up a Polar Operational Environmental Satellites (POES). The framework comprises sets of satellites, which guarantee that each part of the Earth is consistently seen at twice every 12 h from around 800 km elevation. Beginning with the NOAA-15 satellite in 1998, an upgraded form of the Space Environment Monitor (SEM-2) is being flown. The SEM-2 contains two arrangements of instruments that detect the high-energy charged particle counts close to the Earth. It distinguishes and monitors the flux of energetic particles and electrons into the atmosphere and the flux of the particles situated at the height of the satellite. Obara et al. [84] showed the electron flux enhancement due to a geomagnetic storm using an NOAA satellite. Soraas et al. [85] did a similar job using this NOAA satellite data but for the computation of proton flux enhancement related to a geomagnetic storm. Particle bursts due to seismic activity is based on the phenomena of ionospheric-magnetospheric coupling process. The vast description of this coupling process has been given by Walt et al. [86]. Bortnik et al. [87] found the fluctuation in the VLF/LF signal at the same time-frame in which Ref. [86] detected the particle bursts in the radiation belt. Further investigations are conducted by China Seismo-Electromagnetic Satellite (CSES) and ARINA experiments. These examinations show the seismogenic particle bursts [88,89]. Fidani et al. [90] used 0${}^{\circ }$ telescope data from NOAA-15 satellite for 11 years to study the correlation between energetic particle bursts with EQs. They have reported the need for a specific condition to ensure the validity of the decided connections. The relationship is conditioned by several EQs and occasional increments of PBs, which can interfere with each other. Recently, Chakraborty et al. [23] proved significant enhancement of particle count number before two types of the EQ (land and ocean).

In the first 2.5 years of Swarm satellite mission, De Santis et al. [91] studied seismic anomalies before twelve major EQs. They have used a total of 60 days of Swarm satellite data (one month before and one month after the EQ day) for their investigation. For tracking the anomaly, they confined the investigation area by a circle. The EQ epicenter is used as a center for this circle. They have observed fluctuation in the magnetic field and electron density before all these 12 EQs. They also observed that the anomalies reflect a linear relation with the magnitude of the EQ. In the next year, Marchetti et al. [92] reported another case study of EQ anomaly using Swarm satellite data. They have used all the EQs that happened in Central Italy from 2016 to 2017. They also found the anomaly in the Y(East)-component of the magnetic field before the EQs.

In this paper, we adopt a multi-parametric approach to detect the pre-seismic anomalies in the atmosphere and ionosphere for the 2020 Samos EQ. This large EQ, having a magnitude of M = 6.9, occurred in the Aegean Sea, off the coast of the Samos Island (Greece), close to the Greece–Turkey border at 11:51 UTC on 30 October 2020, with a depth of ∼21 km according to the USGS (United States of Geological Survey). In this approach, we have tried to find out seismic anomalies in two different channels of LAIC, namely the acoustic and electromagnetic channels. For the observation purpose, we have used both ground- and space-based instruments. Firstly, for the ground-based observation, we compute the ionospheric TEC anomalies using GPS observation. For the space-based observation, we analyze AGW activity related to the SABER satellite. We also compute the AGW from GPS-TEC to verify the outcomes of SABER. The small-scale function in TEC is detected by the fitting model and wave-like structures are obtained from wavelet analysis. For the ionospheric-magnetospheric coupling process, we examine the perturbation in the variation of particle counts number obtained from the NOAA-15 satellite. From the information of the magnetic field and plasma density from the Swarm satellite, we proceed with a similar methodology as described in De Santis et al. [91] with some necessary modifications to it. In the next section, we present the methodologies we adopt in this manuscript. In Section 3, we present our results, and, finally, in Section 4, we present our conclusions.

2. Methodology

In this manuscript, we investigate the multi-parametric approach for the 2020 Samos EQ, a very strong EQ that took place off the coast of the Northern part of Samos Island (Greece): M = 6.9, epicenter (37.9001${}^{\circ }$ N, 26.8057${}^{\circ }$ E), focal depth = 12 km, time of occurrence 11:51:57 UTC [93,94]. For the ground-based observation, we utilize two GPS-IGS stations (i) DYNG (38.078${}^{\circ }$ N, 23.93${}^{\circ }$ E) in Greece, and (ii) IZMI (38.39${}^{\circ }$ N, 27.082${}^{\circ }$ E) in Turkey, which are close to the EQ epicenter. The locations of the GPS-IGS stations (blue squares), the EQ epicenter (red disk), the EQ preparation zone (EPZ) (blue circle), and the critical zone (CZ) (red circle) are shown in Figure 1. The track of the Swarm satellite is also marked with a black line. The distance of the epicenter from DYNG and IZMI is 251 km and 58 km, respectively. The GPS-RINEX observation and navigation files are taken from the IGS NASA archive (https://cddis.gsfc.nasa.gov/gnss/information/daily, accessed on 22 November 2020). These GPS-RINEX data are fed in a software name “GOPI SEEMALA software” [95] and computed the Vertical Total Electron Content (VTEC) using Equations (1) and (2) [96,97,98].

It is well known that the VTEC can be expressed as [41]:

$$VTEC=\frac{STEC-T_{cal}}{M\left(\alpha \right),}$$

(1)

where STEC is Slant Total Electron Content, $TE{C}_{cal}=b_s+b_R+b_{RX}$, $b_s$ is satellite bias, $b_R$ is receiver bias, and $b_{RX}$ is the receiver inter-channel bias. The $M\left(\alpha \right)$ is obtained by the same methodology as used by Mannuchi et al. [96,97] and Langley et al. [98] as follows:

$$M\left(\alpha \right)=\frac{1}{cos\left(\beta \right)}={(1-\frac{R_{e}cos\left(\alpha \right)}{R_e+h_{min}})}^{-0.5}.$$

(2)

Here, $R_e$ is the radius of the earth and $h_{min}$ is the Ionospheric Pierce Point (IPP) at 350 km [99]. $\alpha$ and $\beta$ are the zenith angle at the receiver site and elevation angle at the IPP, respectively. All of the bias correction and calibration of TEC are made by using the previous methodology in Seemala and Valldares [95] developed and freely distributed by the Institute for Scientific Research, Boston College, MA, USA to compute the STEC and VTEC. It is mandatory to focus on the geomagnetic condition during the investigation of any seismogenic anomalies to eliminate any possible contamination. We gather A_p averages data directly from the World Data Center for Geomagnetism, Kyoto (http://wdc.kugi.kyoto-u.ac.jp/, accessed on 22 November 2020). In Figure 2, we present the variation of the $D_{st}$, $K_p$, $A_p$ 3 h average and daily average, Sudden ionospheric disturbance (SID) and interplanetary magnetic fields (IMF-Bz) during the period from 17 October to 4 November 2020. It is evident that, during that time period, the minimum $D_{st}$ value was −36 nT and the maximum daily sum $K_p$ value was 45 (<50). It is widely published that IMF-Bz and associated variables ($D_{st}$, $K_p$, $A_p$, etc.) are used to indicate the presence of solar-geomagnetic storms. Malik et al. [100] and Ayomide and Emmanuel [101], have analyzed many strong geomagnetic storms with high IMF-Bz values. They confirmed the condition that the value of IMF-Bz within the limit of −10 to +10 nT is considered to be a solar quiet day. They observed a significant enhancement in TEC variation during those strong storms when the interplanetary magnetic fields (IMF-Bz) cross that quiet range. There are also many publications [102,103,104,105,106,107,108,109] where it has been widely found that the moderate IMF-Bz (within the above-mentioned limit) do not have a significant effect on TEC variation.

2.1. Computation of Ionospheric TEC from GPS-IGS Station

To study the anomalies in TEC, we process the median X of the TEC values for the past 15 days of the EQ and the related interquartile range $IQR$. The upper bound ($UB$) and the lower bound ($LB$) at a specific time (UT) is given by

$$UB=X+IQR,$$

(3)

$$LB=X-IQR.$$

(4)

We use a technique similar to that used by [28] to calculate the anomalies in VTEC. Assuming a conventional dispersion with a mean ($\mu$) and standard deviation ($\sigma$) for the VTEC, the usual values of X and $IQR$ are $\mu$ and 1.34$\sigma$, respectively [110]. The VTEC crosses either the lower or upper bound with an 80–85% confidence level that is chosen as an anomaly. After calculation of the upper and lower bound, we determine the anomaly by taking the enhancement and decrements above and below the upper and lower bound, respectively [52]:

$$anomaly=VTEC\ge UB$$

(5)

$$anomaly=VTEC\le LB$$

(6)

2.2. Computation of AGWs from the SABER Temperature Profile

SABER is an instrument of the TIMED satellite that was launched in December 2001 at an orbital height of 625 km. The period of the satellite is 1.7 h with a tendency at 74.1${}^{\circ }$ [111]. SABER mentioned its first observable fact in January 2002. It records the temperature with a height inclusion of 20–100 km by utilizing the wavelength range from 1.27 to 17 $\mathsf{\mu }$m. From toward the north view, the latitudinal scope of SABER is between 50${}^{\circ }$ S–82${}^{\circ }$ N and from toward the south view, it is 50${}^{\circ }$ N–82${}^{\circ }$ S. The perception seeing method of SABER changes in each 60-day time frame. Remsberg et al. [112] have examined the technique for extraction of temperature from the SABER satellite.

Preusse et al. [113], Preusse et al. [114], and Fetzer et al. [115] have effectively built up the strategy to compute the AGW. The mechanism to get the AGW from the temperature profile has been acquired by a few techniques in the previous few decades [21,22,71,73,116,117,118,119,120,121,122]. In the following, we summarize the methodology applied for the extraction of the AGW from the SABER temperature profile. We gather the altitude variation of temperature profile for the elevation range 20 to 100 km for the region of our investigation (around the EQ epicenter) from the SABER archive (http://saber.gats-inc.com/, accessed on 22 November 2020) and take the logarithm of the acquired individual temperature profiles. At that point, a third-order polynomial is fitted on the logarithm temperature profile and, to get the residual temperature, we subtract the fitted profile from the original one. As AGW has a wavelength longer than 4 km, so to eliminate the other small-scale waves, a 4 km boxcar filter is applied to the residuals of the individual profiles. After such filtration, the filtered data are added back to the fitted profile, and this gives the final profile. An antilogarithm of the final profiles is known as least square fit (LSF), which is utilized to acquire the daily zonal mean temperature and other zonal wave components from 1–5. The obtained background temperature ($T_0$) is obtained from the summation of all wave components 0–5. The background temperature profiles are subtracted from the original temperature profile to get the perturbation temperature ($T^{\prime }$). The obtained background temperature profiles are put in Equation (8) to get the Brunt Vaisala frequency (N) for a particular profile:

$$N^2=\frac{g}{T_0}(\frac{\partial T_0}{\partial z}+\frac{g}{c_p}),$$

(7)

where z is the altitude, and $c_p$ is the specific heat at constant pressure.

The potential energy ($E_p$) associated with the AGW can be estimated for individual temperature profiles by putting the perturbed and background temperature in Equation (6):

$$E_P=\frac{1}{2}{(\frac{g}{N})}^2{(\frac{T^{\prime }}{T_0})}^2$$

(8)

where g is the acceleration due to gravity, and N is the Brunt Vaisala Frequency.

This technique is now utilized in our past attempts to contemplate the activity of AGW during an EQ that happened on 3 January 2016, in Imphal, India [73] and the investigation of AGW action during numerous EQs. The strategy is initially utilized by [119] to remove the GW from the SOFIE temperature profile to examine the activity of AGW during different Sudden Stratospheric Warming (SSW) events. To compute the AGW excitation during the EQ, we choose a period from 17 October to 4 November 2020. The AGW is computed for a region of 30${}^{\circ }$–50${}^{\circ }$ N (latitude range) and 10${}^{\circ }$–40${}^{\circ }$ E (longitude range), which mostly covers the EQ epicenter and the EPZ (∼1023 km) and CZ (∼200 km).

2.3. AGW Detection from GPS-TEC

Alternatively, we use the GPS-TEC information to get the AGW excitation that is a TID of a period between 60 to 120 min. Fundamentally, the perturbations made by the TIDs are small-scale fluctuations that are essentially produced because of the EQ and meteorological events. This small-scale fluctuation can be captured by GPS signals during proliferation. Even at a short distance, the occurrence of TID often leads to large gradients in the TEC. To detect this, we use a fitting method known as “Savitzky–Golay filtering” (sgolayfilt) [123,124] (Equation (9)). In Savitzky et al. [123], they showed that the weighting coefficients for a smoothing operation may be determined by a set of integers $(C_{-n},C_{-\left(n1\right)}...,C_{n-1},C_n)$. The utilization of these weighting coefficients, known as convolution integers, is precisely the same as fitting the information to a polynomial, as explained and is computationally more efficient and faster. The Savitzky–Golay technique therefore provides the smoothed data point ${\left(x_i\right)}_s$ using the following equation:

$${\left(x_i\right)}_s=\frac{{\sum }_{i=-n}^{i=n}{C}_i{x}_{i+1}}{{\sum }_{i=-n}^{i=n}{C}_i}$$

(9)

The deviation of $dVTEC$ is obtained by subtracting the modeled profile from the original data. To obtained the modeled $VTE{C}_f$ profile, we utilize the sgolayfilt with polynomial order 5 and 90 min time window in MATLAB. We use the filtering in a large time window to remove the additional noise in the $dVTEC$. The small-scale fluctuations are computed by taking the difference between the observed and fitted values as shown in Equation (10) [58]:

$$dVTEC=VTEC-VTE{C}_f$$

(10)

After getting the small-scale fluctuations ($dVTEC$), a spectral analysis of $dVTEC$ is performed to examine the possible wave-like structures associated with it. We do a wavelet scalogram analysis of $dVTEC$ by using the complex Morlet continuous wavelet through MATLAB [125]. Given the computation of the scalogram, we investigate in more detail wave movement in the time frame of 15 min to 1 h (i.e., MSTIDs) and of 30 min to 3 h (i.e., LSTIDs). The maximum power of the spectrum (MPS) represents proxies for the typical spectral amplitude for a given measurement path in the individual range. The Wave like structures within the Cone of Influence (COI) are the signatures of the AGWs [58].

2.4. Computation Process of Energetic Particle Bursts

It is notable to us that a large EQ greatly affects radiation belt energetic particle counts and, for this reason, we have assembled NOAA satellite information. NOAA-15 satellite gives processed and raw data information, and we choose these datasets for a period of 26 days starting from 15 October 2020 to 9 November 2020. For the particle counts calculation, our initial step is to make another dataset on an everyday schedule. Each dataset contains time in milliseconds, latitude in degree, longitude in degree, MEPED (Medium Energy Proton Electron Detector) electron channel information (electron count rates for example CR’s), IGRF magnetic field (B), MEPED telescope pitch angle ($\alpha$), and L values (McIlwain L-parameter). We average out these datasets every 8 s. To eliminate this cumulative sum in the newly-made datasets, we make a difference in the energy values and growing new energy channels according to our further computation. These new energy channels are 30–100 keV, 100–300 keV, and <300 keV; likewise, we exclude the energy values that contain zero. We make a three-dimensional matrix with L, $\alpha$, and B values where the binning of L value is from 0.9 to 2.2 with a bin width of 0.1. The range of the pitch angle (0${}^{\circ }$ to 180${}^{\circ }$) is divided with 12 stages of 15${}^{\circ }$ each. Geomagnetic field B is separated into nine ranges as follows: 16.0–17.5 $\mathsf{\mu }$T, 17.5–19.0 $\mathsf{\mu }$T, 19.0–20.5 $\mathsf{\mu }$T, 20.5–22.0 $\mathsf{\mu }$T, 22.0–25.0 $\mathsf{\mu }$T, 25.0–29.0 $\mathsf{\mu }$T, 29.0–33.0 $\mathsf{\mu }$T, 33.0–37.0 $\mathsf{\mu }$T, and 37.0–41.0 $\mathsf{\mu }$T and, by the utilization of this matrix, we need to count the high energetic particles included in each shell distinguished by the NOAA-15 satellite. Now, our next step of this computation is to make a condition for which the count rates are considered as particle bursts (PB’s). Based on the previous theory, the 8 s CR’s are consistent with Poisson distribution [82]. Thus, to determine the non-Poissonian fluctuation with 99% probability, the defined condition is 4$\sigma$ value. The energetic particle numbers which exceed this value are considered PB’s. Now, the seismic-induced PBs are selected by putting up the condition $\left|\mathsf{\Delta }L\right|$≤ 0.1, where $\mathsf{\Delta }$L is characterized by the contrast between the L values related to EQ and particle bursts.

2.5. Analysis of Swarm Satellite Data

The fifth mission in ESA’s fleet of Earth Explorers is Swarm that provides magnetic field and plasma information by its three satellites. These are satellite ALPHA (SAT-A), satellite BRAVO (SAT-B), and satellite CHARLIE (SAT-C). It was launched on 22 November 2013, and the final constellation was achieved on 17 April 2014. The inclination angle is 87.35${}^{\circ }$ for SAT-A and SAT-C at an altitude of 462 km (initial altitude), and the SAT-B is inclined at 87.75${}^{\circ }$ and having a higher altitude (511 km) than the other two satellites. SAT-A and SAT-C are flying with a fixed longitude separation of 1.4${}^{\circ }$. Swarm consists of an absolute scalar magnetometer (ASM) which can measure the strength of the magnetic field and provides the scalar value of the magnetic field and an electric field instrument (EFI) having two components i.e., Langmuir Probe (LP) and Thermal Ion Imager (TII). Electron density, electron temperature, and spacecraft potential data are provided by LP, and TII provides ion drift and ion velocity in higher resolution. Now, the third instrument in Swarm is the vector field magnetometer (VFM), which gives a magnetic field in vector form. This instrument also provides altitude data. Besides these three main instruments, Swarm has a star tracker, GPS receiver, and an accelerometer.

Magnetic Field and Plasma Data Structure in Swarm Satellite

VFM provides the magnetic field data that consist of magnitude and direction because it is a fluxgate magnetometer with a compact spherical coil sensor, and it gives data with two types of resolution. One is a higher resolution having a sampling rate of 50 Hz and another one is a lower resolution with a sampling rate of 1 Hz. In a reference frame of North, East, and Center (NEC), the geomagnetic field values are provided in a spherical coordinate. ASM is based on the electron spin resonance (ESR) theory using the Zeeman effect, and it is used to calibrate the VFM data. Swarm data are available in ESA earth online and are freely accessible (http://Swarm-diss.eo.esa.int, and ftp://Swarm-diss.eo.esa.int, accessed on 22 November 2020). In our study, we use low rate (1 Hz) VFM data, and we use time (in UTC), latitude (in degree), longitude (in degree), X-component, Y-component, Z-component data of magnetic field strength, scalar intensity (F), and also the plasma density, local temperature of plasma, and the potential energy of spacecraft. We mainly use electron density as a vital parameter to detect the anomaly due to EQs and want to verify the possible variation in this parameter related to the 2020 Samos EQ. For computation purposes, we follow three algorithms by which we can easily determine the unusual variation in the magnetic field component and in the electron density which are corresponding to our study.

• I. MASS Algorithm:

Anomalies in the magnetic field component are searched in each track of the Swarm satellite by MASS (magnetic Swarm anomaly detection by spline analysis). It follows a few steps such as:

Firstly, we summarize our conventional data that is being used in the MASS algorithm from CDF data. After extracting these CDF datasets, the extracted data file consists of additional information of SAT-A, SAT-B, and SAT-C, semi-orbital track number, data quality flag, time (in UTC), and local time (LT). From the extracted data, we create a file having seven columns (time (UTC), latitude (in degree), longitude (in degree), scalar form of the magnetic field intensity, X-component of the magnetic field, Y-component of magnetic fields, and Z-component of the magnetic field. For the Charlie satellite, no ASM (scalar) data are available from 5 November 2014 (7:37 p.m. UTC) due to a technical problem. Thus, we use only VFM data for SAT-C. Next, we transform all the geographic latitude into geomagnetic latitude using IGRF-12 (https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html, accessed on 22 November 2020).

It is well known that geomagnetic storms play a major role on ionospheric irregularities. Thus, for the accuracy of the anomaly related to our considering EQ, we exempt the geomagnetic effect on the entire calculation. To eliminate the geomagnetic activity from this computation process, we analyze the geomagnetic index Dst and A_p values for each time of the different satellite tracks (every one hour of Dst and every 3 h of Ap). The analyzed area is already mentioned in SABER methodology. Now, we develop different codes in MATLAB following these steps:

We plot each satellite track from 0 h to 24 h. We consider only those tracks that pass very close to the EQ epicenter to detect the anomaly related to this EQ. As the epicenter of the 2020 Samos EQ is 37.918${}^{\circ }$ N and 26.79${}^{\circ }$ E, we chop the Swarm data in the latitude range 0${}^{\circ }$ to 55${}^{\circ }$ N and longitude range 0${}^{\circ }$ to 55${}^{\circ }$ N and finally create the individual dataset for each hour from 0 to 24. We take the first derivative (difference between two consecutive values) of the magnetic field components in each track to acquire more information about the data and to remove the higher degree values. Now, we use the cubic spline method as the best fitting technique to remove the remaining long trend along the selected track. We apply the fit to the X, Y, Z component of the magnetic field and F. Finally, we compute the residue corresponding to the best fit and identified the anomaly related to the seismic event.

• II. NeLOG Algorithm:

Electric field instrument data provide the electron density with a sampling rate of 2 Hz. By the use of the NeLOG program, we successfully identify the anomalies in the electron density on a logarithmic scale. Regarding 2020 Samos EQ, we analyze the electron density, electron temperature, and potential energy data from one month before to one month after the EQ day. In the NeLOG program, we select the ±15${}^{\circ }$ around the epicenter and time is taken from 12:00 p.m. UT to 1:00 p.m. UT. This algorithm provides the plot of decimal logarithmic electron density, and this plot is compared with the geographic representation of the EPZ, CZ, with EQ epicenter and the selected track. We only use electron density values as it is less affected by the instrumental error than the other plasma parameters. Now, we fit the plot with a 10-degree polynomial as it is the best fit for this variation of logarithmic electron density. Finally, the root means square is calculated within the above-mentioned area. To eliminate the fitting edge error, we cut the first 5${}^{\circ }$ from the fitted curve. In our case, we choose track number 15 as it contains all the anomalous behavior related to the EQ. We only consider SAT-C because SAT-A and SAT-B do not provide electron density information from 24 to 26 October 2020.

• III. NeSTAD Algorithm:

The Ne single track anomaly detection (NeSTAD) algorithm is used to detect the anomalous behavior of electron density in a given interval of time. For our case, this interval is between 12:00 p.m. UTC to 1:00 p.m. UTC on 29 October 2020 (one day before the 2020 Samos EQ). This program runs with the input data of Langmuir Probe of Swarm satellite (Swarm EFIx-PL and Plasma Preliminary). We first select the region for examining the geographical range as mentioned in the MASS algorithm and similar binning for the latitude–longitude in the above-mentioned time interval. To remove the large-scale electron density gradient in other regions, we select the geographic region for our analysis to determine the anomaly in the electron density data in a single track. Now, for every single track in the geographic range and the time range, the following parameter is calculated:

$$\mathsf{\Delta }{\left(\frac{\mathsf{\Delta }{N}_e}{N_e}\right)}_{ti}=\left(\frac{N_e\left(t_i\right)}{N_e\left(t_{i+1}\right)}\right)-\left(\frac{N_e\left(t_{i+1}\right)}{N_e\left(t_{i+2}\right)}\right).$$

(11)

An outlier is simply an unusually large or small value among the rest. In statistical analyses, they may generate issues. We computed the outliers in the anomalous track with respect to the normal tracks. For this computation, we have used a ±15 day period from the day of the EQ. It is necessary to utilize interquartile ranges (IQR) to set the outliers. This is the difference between the 3rd and 1st quartiles (IQR = Q3 − Q1). Outliers are defined in this context as the distribution $\mathsf{\Delta }\left(\frac{\mathsf{\Delta }{N}_e}{N_e}\right)$ are either below (Q1 − k· IQR) or above (Q3 + k·IQR). The following quantities (referred to as fences) are required to identify extreme values in the distribution’s tails for the values k = 1.5 or 3: lower inner fence is Q1 − 1.5· IQR, upper inner fence is Q3 + 1.5· IQR, lower outer fence is Q1 − 3· IQR, and upper outer fence is Q3 + 3· IQR. A mild outlier is defined as a point beyond an inner fence on each side. An extreme outlier is a point that lies beyond the outer fence [126] (NIST/SEMATECH e-Handbook of Statistical Methods, http://www.itl.nist.gov/div898/handbook/, accessed on 22 November 2020).

Finally, $\mathsf{\Delta }{\left(\frac{\mathsf{\Delta }{N}_e}{N_e}\right)}_{ti}$ is calculated for each track of Swarm satellite. The quantity $\mathsf{\Delta }{\left(\frac{\mathsf{\Delta }{N}_e}{N_e}\right)}_{ti}$ is able to identify steep variations in the electron density time profile, behaving like a high-pass filter on the electron density data. The track anomaly parameters are derived by

$$R=1-max\left(\mathsf{\Delta }\_filt\right))/max\left(\mathsf{\Delta }\right),$$

(12)

where $\sigma$ is the standard deviation of $\mathsf{\Delta }\left(\frac{\mathsf{\Delta }{N}_e}{N_e}\right)$ and % = percentage of outliers recognized in the track, $\mathsf{\Delta }\_filt$ denotes the strength of the outliers after the filtration process, and $\mathsf{\Delta }$ denotes the strength of the outliers before the filtration process.

3. Results

3.1. Ionospheric Perturbation Observed from GPS-TEC

The diurnal variation of VTEC from 17 October to 4 November 2020, with the upper and lower bound are shown in the upper panel of Figure 3 for the DYNG station. The maximum enhancement in VTEC is observed on 29 October. In addition, on 22 and 24 October, similar unusualness is observed. The lower panel quantifies the fluctuations in TEC. It is found that the anomaly in TEC is 2.5, 2.9, and 4.5 TECU on 22, 24, and 29 October, respectively. It is evident that the pre-seismic anomalies started eight to nine days before the EQ and becomes the maximum just a day before the EQ. Figure 4 is the same as Figure 3 for the IZMI station. For the IZMI station, the daily TEC variation crosses the upper bound on 21 and 29 October with an amount of anomalies of around 2.5 and 3 TECU, respectively. Therefore, for both of the stations, the VTEC is only the maximum on just the day before the EQ. Even though IZMI is closer to the epicenter, maximum changes in TEC are observed in the DYNG station. It can be concluded that our anomalous TEC variation during the EQ was not influenced by these low IMF-Bz values (22 October 2020: 8 nT, 24 October 2020: 11 nT, 29 October 2020: 6 nT). Though, on 24 October 2020, the IMF Bz maximum value is slightly greater than 10, it could not have a significant impact on TEC variation that can contaminate the anomalies in TEC from the seismogenic cause.

3.2. AGW Anomalies Observed SABER Satellite

After the computation of $E_p$ (See Section 2.2), a nine-dimensional matrix is obtained for individual temperature profiles. The nine-dimensional matrix contains latitude, longitude, date, or day of the year in UT, altitude, original SABER temperature profile, reconstructed fitted temperature profile, perturbation temperature, Brunt Vaisala frequency, and AGW associated potential energy ($E_p$).

The time–altitude variation of $E_p$ with the associated AGW from 17 October to 4 November is shown in Figure 5. The altitude ranges from 30 to 50 km. The $E_p$ is presented as a color bar that ranges from 0 to 10 J/kg. This spatio-temporal profile indicates a significant amount of $E_p$ during 23 October to 25 October at around 46 km to 48 km altitude. The EQ day is marked with a black dashed line, and it is evident that the AGW activity is significantly enhanced around six to eight days before the EQ.

Figure 6 describes the spatial distribution of $E_p$ for the same time period. Here, the intensification of $E_p$ is projected on a two-dimensional map for which the latitude and longitude range from 30${}^{\circ }$ N to 45${}^{\circ }$ N and 10${}^{\circ }$ E to 40${}^{\circ }$ E with the EQ epicenter at approximately in the center (magenta diamond). We take the altitude of 47 km and present the $E_p$ values. Our selection of altitude follows the maximum enhancement in $E_p$. It is evident from Figure 6 that maximum AGW activity is observed on 24 October 2020, near the EQ epicenter. The patch of AGW stays at the northeast of the epicenter and the maximum enhancement is found within the 500-km radius from the epicenter. On 25 October, a similar but much smaller patch is observed in an opposite direction of the epicenter.

It is obvious that there is a strong signature of the presence of AGW as computed from the SABER satellite outcomes. There is no such secondary strong peak of AGW in both temporal and spatial variation of AGW, and thus it can be generated due to the studied EQ.

3.3. AGW Anomalies from GPS-TEC

To validate the outcomes of AGW from SABER observation, we present the indirect method of AGW computation as retrieved from the GPS-TEC findings. For TEC fluctuations, the normal unperturbed ionospheric condition follows the range of −0.25 ≥ dVTEC ≤ 0.25. Any perturbation for which the estimation of dVTEC or the small scale fluctuations in dVTEC value crosses this envelope is considered as a seismogenic anomaly. It is evident from Figure 7 that the estimation of dVTEC satisfies the above-mentioned condition for 19, 21, 27, and 29 October for the DYNG station. The maximum fluctuation dVTEC is seen on 19 November 2020.

To examine the presence of possible wave-like structures in the similar temporal domain, the corresponding scalogram is illustrated in Figure 8. From the scalogram, a wave-like structure of period 65 to 80 min appears on 19 October 2020 inside the COI. A similar enhancement in terms of wavelet periods ranging from 65 to 90 min and 65 to 100 min is also visible on 21 and 29 October 2020 with less intensification. It additionally showed a wave-like structure on 21 and 29 October 2020 with the time of 65–90 min and 65–110 min, respectively. The maximum intensification of 19 October is plotted separately in Figure 9 for better representation. The top, middle, and bottom panels show the observed and fitted TEC profile, dVTEC, and the scalogram, respectively. The scalogram shows significant enhancement of wave-like structure having the period of gravity waves within the COI.

Figure 10 and Figure 11 present the similar computation of dVTEC and corresponding scalograms for IZMI station, respectively. In contrast to the DYNG station, there are significant fluctuations in dVTEC from 19 to 22 October, 25 to 26 October, and 29 October. The corresponding wavelets corroborate the same where an intense wavelike structure is observed for those similar days. Just like IZMI, the intensification is a maximum on 19 November with a period of 60 to 80 min. A sharp difference is observed in the scalogram IZMI where the intense path of AGW has more temporal spread over the day. As similar to Figure 9, the overall scenario of 19 November is shown in Figure 12.

3.4. Anomalous Energetic Particle Burst (PB)

Figure 13 shows the daily average count rates (CRs) for 30 October 2020 (day of the EQ). Each 3D shell reflects the CRs for which the satellite passes at least 20 times. The x-axis indicates the L values, the y-axis shows the pitch angle, and the color bar indicates the number of times the satellite passes in each 3D shell. The satellites pass ≥20 times for a cell for which the magnetic field value exceeds 22.0 $\mathsf{\mu }$T. As mentioned in the previous section, we choose those regions for further analysis, for which this magnetic field condition is satisfied. In our investigation, we have barred the South Atlantic Anomaly (SAA) and the external VAB by picking B > 22.0 $\mathsf{\mu }$T also, L < 2.2 individually. By this procedure, we eliminate the SAA part from the entire analysis of particle bursts computation, as the South Atlantic region always shows a sharp gradient of particle counts number.

CRs that exceed four times the standard deviation, according to the Poisson distribution, are detected as non-Poisson fluctuations with a probability of 99%, and these CRs are called particle burst events [82,83]. To eliminate the normal background statistical fluctuation, we set a constraint to eliminate the same by computing the $\sigma$ values of the particle counts profile. In Figure 14, we present the 8 second average of particle counts and the computed 4$\sigma$ level. The particles that crossed this level of 4$\sigma$ are treated as the particle bursts, and the remainder of the number of particles are viewed as the normal statistical fluctuation. Here, the x-axis demonstrates the electron included in 8 s average and, along with the y-axis, it shows the daily energetic particle count numbers. The blue dashed line demonstrates the 4$\sigma$ level. For this EQ, the estimation of this 4$\sigma$ is 3.32.

Energetic particles are getting perturbed in the radiation belt due to many other sources. Thus, the elimination of these sources is an essential task to take care of additional examination for the true detection of EQ-associated particle bursts. Thus, to take out the other sources which assume a vital part to enhance the quantity of the high energy particle count number, we examine the A${}_p$ index values which genuinely distinguish the presence of solar activity. We utilize the geomagnetic record A${}_p$ information from the World Data Centre for Geomagnetism, Kyoto. From the American Association of Variable Star Observers (AAVSO) (https://www.aavso.org/sid-database, accessed on 06 July 2021) database, solar quiet days were characterized as those with daily averages of A${}_p$ < 16, A${}_p$ < 25 in each of the three hour intervals throughout the day and SID = 0. If any of these conditions were violated on any given day, we designated that day as a solar active day (A${}_p$ < 16, A${}_p$ > 25, SID = 0 or A${}_p$ < 16, A${}_p$ < 25, SID = 1) [82,83,127]. Based on this selection, we separate the contaminated and non-contaminated PBs and effectively eliminate the geomagnetic impact on the high energetic particle count number. We have found the solar active days are 290, 292, 297, 298, 301, 303, 306, 309, 310, 311, 312, and 313 (red histogram shown in Figure 15 upper panel). The rest of the days are solar quiet days (black histogram shown in Figure 15 upper panel). The monthly quiet day average value of PB is 30.4555, denoted by a horizontal dashed line. As mentioned earlier (See Section 2.4), by putting $\left|\mathsf{\Delta }L\right|$≤ 0.1 condition, we successfully detected the seismic-induced PBs for the studied EQ. The lower panel of Figure 15 shows that the number of PB is zero on the day of the EQ (30 October 2020). We have observed that the EQ induced PB on the day of years 294, 300, and 305. There are a significant number of PBs on ten, four days before the EQ. It is obvious that, even though there are few PBs found in this entire observation period of 26 days, not all PBs are initiated by this EQ. The lower and upper panel of Figure 15 show that mismatch. The reason for such decrement (in lower panel) in the number of PB is the removal of PBs in our computation that originates from various L-shells. We only consider those PBs which are produced from the L${}_{eq}$ (EQ-related L-shells). It is also evident that, for a geomagnetically quiet condition, the EQ-induced PBs (See Figure 15 lower panel) have a lower value than those number of PBs for which we do not impose the $\mathsf{\Delta }$L condition although they are non contaminated (see Figure 15 upper panel).

3.5. Outcomes from the Swarm Satellite

From the MASS algorithm, Figure 16 shows that the anomaly is mainly observed in the X component of the magnetic field (first panel) in SAT-C. This disturbance occurred one day before (29 October 2020) the EQ at the geomagnetic latitude range 10${}^{\circ }$ N to 12${}^{\circ }$ N, while the source (EQ epicenter) is at the geomagnetic latitude 32.20${}^{\circ }$ N. No such anomaly is observed in the Y and Z components of the magnetic field (second and third panels). As SAT-C does not provide the scalar component values of the magnetic field, it remains constant (fourth panel). The anomalous day is geomagnetically quiet having A${}_p$ = 27 nT, K${}_p$ = 4, and D${}_{st}$ = −28 nT. This anomalous track is detected in the afternoon time 1:00 p.m. UTC (Local Time = 15 EET), which is a strong and unexpected behavior in the magnetic field lines, as it was geomagnetically quiet.

The anomalies in the electron density profile from the NeLog algorithm are shown in Figure 17. Figure 17 shows the comparative results of the MASS algorithm (first to the fourth panel) along with the result of the NeLog algorithm (fifth panel) for the same day and same time as mentioned above. We observe small fluctuations around 8${}^{\circ }$ N to 11${}^{\circ }$ N geomagnetic latitude in the electron density profile (fifth panel). In Figure 18, we present Log${}_{10}$ (N${}_e$) (first panel), time derivative of electron density (second panel), and the residual value of the electron density (third panel). The electron temperature (black line) and the corresponding potential energy (red line) are presented in the fourth panel, respectively. The parameter d$N_e$/dt (time derivative of electron density) confirms an anomaly in the latitude range 8${}^{\circ }$ N to 11${}^{\circ }$ N. A similar kind of fluctuation is observed in electron temperature (T${}_e$) in the same latitude range Figure 18.

The observation through NeSTAD analysis is shown in Figure 19 for the SAT-C. Figure 19 represents the electron density variation (left panel) with geographical latitude range from 5${}^{\circ }$ N to 30${}^{\circ }$ N at 13:00 UTC. In the right panel, the red curve represents the data of $\mathsf{\Delta }\left(\frac{\mathsf{\Delta }{N}_e}{N_e}\right)$ before the outlier analysis and the black curve denotes the variation after the outlier analysis. The value of the “track anomaly parameters” are $R=0.957$, $\sigma =0.0085$, and $\%=0.0423$, respectively.

Figure 20 shows a comparison of dates of abnormal activity in the stratosphere, ionosphere, and magnetosphere related to 2020 Samos EQ (

Table 1. Table for acronyms.

Full Form	Acronym
Earthquake	EQ
Lithosphere–Atmosphere–Ionosphere Coupling	LAIC
Atmospheric Gravity Wave	AGW
Total Electron Content	TEC
Low Earth Orbit	LEO
National Oceanic and Atmospheric Administration	NOAA
Sounding of the Atmosphere using Broadband Emission Radiometry	SABER
Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics	TIMED
Count Rate	CR
South Atlantic Anomaly	SAA
Particle Bursts	PB
Very-Low-Frequency	VLF
Low-Frequency	LF
Ultra-Low-Frequency	ULF
Detection of Electromagnetic Emissions Transmitted from Earthquake Regions	DEMETER
CHAllenging Minisatellite Payload	CHAMP
Global Positioning System	GPS
Outgoing Longwave Infrared Radiation	OLR
Atmospheric Chemical Potential	ACP
Frequency Magnitude Relationship	FMR
Global ionosphere map	GIM
Traveling ionospheric disturbances	TIDs
Global Assimilative Ionospheric Model	GAIM
Gravity Wave	GW
Middle and upper environment	MU
Mesosphere–Stratosphere–Troposphere	MST
Van Allen Radiation Belt	VAB
National Aeronautics and Space Administration	NASA
Polar Operational Environmental Satellites	POES
Space Environment Monitor	SEM
China Seismo-Electromagnetic Satellite	CSES
EQ preparation zone	EPZ
Critical Zone	CZ
Sudden Ionospheric Disturbance	SID
Interplanetary Magnetic Fields	IMF-Bz
Vertical Total Electron Content	VTEC
Interquartile Range	IQR
Upper Bound	UB
Lower Bound	LB
Least Square Fit	LSF
Savitzky–Golay filtering	sgolayfilt
Maximum Power of the Spectrum	MPS
Cone of Influence	COI
Medium Energy Proton Electron Detector	MEPED
International Geomagnetic Reference Field	IGRF
Langmuir Probe	LP
Thermal Ion Imager	TII
Absolute Scalar Magnetometer	ASM
Vector Field Magnetometer	VFM
Thermal Ion Imager	TII
North, East, Center	NEC
Electron Spin Resonance	ESR
Magnetic Swarm Anomaly Detection by Spline Analysis	MASS
Ne Single Track Anomaly Detection	NeSTAD
Surface Latent Heat Flux	SLHF
Relative Humidity	RH
Sudden Stratospheric Warming	SSW
Solar Occultation for Ice Experiment	SOFIE
Receiver Independent Exchange Format	RINEX
Magnetic Swarm Anomaly Detection by Spline Analysis	MASS
American Association of Variable 645 Star Observers	AAVSO

) as indicated by the multiparametric analysis technique. This figure indicates that AGW activity in the stratosphere was reached around six days prior to 2020 Samos EQ. However, the TEC anomalies that were observed on 8, 6, and 1 day before the EQ day occurred over a longer period of time, with the Swarm magnetic field and plasma density being observed over the course of 1 day before the EQ in ionosphere. This causes the PB, as the resulting EQ, to be observed for 10 and 4 days before the EQ day. On days 10 and 4 before the day of the EQ, the magnetospheric perturbation (here PB) caused by this event has been noticed.

4. Discussion

This work presents a multi-parametric approach to study the pre-seismic anomaly during and before the 2020 Samos EQ that took place on 30 October 2020 in Greece. To achieve our goal, we use a group of ground and space-based techniques and observe stratospheric, ionospheric, and magnetospheric parameters.

• At first, regarding ground-based observation, we use ionospheric GPS-TEC information from two IGS station DYNG (Greece) and IZMI (Turkey) which are close to the EQ epicenter. We compute the diurnal TEC variation and, to detect the pre-seismic anomaly in it, we use the method of statistical upper and lower bounds. The pre-seismic enhancement starts around 8–9 days before the EQ, and the maximum enhancement occurred one day before the mainshock for both DYNG and IZMI stations. The maximum anomaly in TEC is found of 4.5 TECU for DYNG station, which is comparatively farther from the epicenter, while the maximum change is found to be 2.5 TECU for IZMI.
• For computing the AGW associated with the EQ, we use both direct and indirect methods. In the direct method, the AGW activity is observed through the space-based satellite SABER/TIMED. We computed the potential energy $E_p$ associated with the AGW from the atmospheric temperature profile as recorded from SABER. A significant enhancement in $E_p$ associated with AGW is observed 6–8 days (23 to 25 October 2020) before the EQ. The enhancement of $E_P$ is found to be at an altitude range of around 46–48 km. From the spatial variation of AGW, we observe the maximum enhancement in $E_p$ on 24 October 2020, at 47 km altitude with a radius of 500 km in the northeast direction of the EQ epicenter.
• To validate the outcomes of SABER, we use another indirect method, where wave-like structures are investigated in the small scale fluctuations from GPS-TEC using a filtration method. The wave-like structures of periodicity 65–110 min are obtained from the wavelet analysis of small-scale fluctuations for both of the stations. We observe the most intense wave-like structure on 19 October 2020, for both stations. The AGW enhancement in the wavelet spectrum is much concentrated for the DYNG station, whereas, for the IZMI station, it is scattered around a period.
  In space-based observation, we use two satellites, namely NOAA-15 and Swarm, to investigate ionospheric and magnetospheric irregularities associated with the EQ.
• Based on the NOAA-15 satellite particle database, we computed radiation belt energetic particle counts associated with the EQ. By eliminating SAA and considering geomagnetic quiet conditions, we present the number of particle bursts. We observe a significant number of particle counts on 10 and 4 days before the EQ.
• We examine the anomalies in the magnetic field, electron temperature, and electron density profile using Swarm satellite magnetic field and plasma information. In the computation, we use MASS, NeLOG, and NeSTAD algorithms (demonstrated in the methodology section). We observe the anomalous behavior in the X component of the magnetic field by using the MASS algorithm. The anomalous track number is 15 (SAT-C), and the fluctuation in the magnetic field is observed one day before the EQ. This fluctuation was observed at the afternoon period (12:00 p.m.–1:00 p.m. UT) and around −15${}^{\circ }$ latitude from the epicenter. We also observe the anomaly in the time derivative of electron density and electron temperature, around the same latitude and at the same period as derived from the NeLOG algorithm. In addition, by using the NeSTAD algorithm, we observed a similar anomaly in the strength of the outliers. For this case, we detected the mild outliers having a k value of 1.5, associated with this EQ. We recognize 0.0423% of outliers in the anomalous tracks.

Our investigation based on this multi-parametric approach re-established some of the major facts of the LAIC mechanism. First of all, the pre-seismic processes (anomalies in diurnal TEC, enhancement in AGW activity, generation of the EQ-induced PBs, anomalies in the magnetic field, and electron density) are the key ingredients of the three channels of LAIC having completely different characteristics. Although all of the parameters are found to be detectable before the EQ, they have shown similar kinds of differences following the precursory time profile. For instance, the magnetic field anomaly, rate of change of electron density, electron temperature, electron potential, and GPS-TEC showed an anomaly for a short period before the main-shock. These parameters become maximum just one day before EQ day studied. For the other parameters like particle precipitations, the pre-seismic irregularities have an intermediate time distance from the EQ day. Most importantly, as the time period $\pm 15$ days around the EQ is geomagnetically quiet, for the obvious reason, the anomalies observed in all the chosen parameters are possible due to the EQ. We also examined some other thermal parameters such as the surface latent heat flux (SLHF), relative humidity (RH), and outgoing long-wave radiation (OLR), which are usually found anomalous but could not identify any significant indication for this EQ.

5. Conclusions

This manuscript deals with pre-seismic anomalies during the 2020 Samos EQ as observed from some well-known parameters as observed and computed from ground and satellite sources. We examine two major channels viz. (a) acoustic and (b) electromagnetic of the LAIC mechanism. We use TEC, AGW, energetic particle burst in radiation belt, magnetic field, electron density, and electron temperature for our study. In terms of analytical and observational points of view, these parameters are found to be successful for a convincing pre-seismic signature; however, there is a significant difference in a precursory time frame among them. The overall parameters show the pre-seismic anomalies from ten to one day before the EQ. We follow the EPZ concept, and all the parameters are taken within that zone. It is evident from Figure 2, however, that there are few examples of quasi-active days in between the solar quiet days where the A${}_p$ values marginally touch the upper threshold values. We eliminate the possible contamination by using the proper choice of such variables carefully to get only the seismogenic irregularities in a solar quiet condition. This will eventually minimize the possibilities of any such contamination in the seismogenic impression.

It is well established that the wide range of pre-seismic processes perturbs the atmosphere and ionosphere with different spatio-temporal ranges. Thus, it is expected that the parameters of each domain of the LAIC (beneath the surface of the earth, surface of the earth, troposphere, stratosphere, ionosphere, and magnetosphere) will excite and show intensification from their normal value with a different time range. In addition, it is extremely important to understand the preparation mechanism of an earthquake in terms of the physical processes beneath the earth that deal with the generation of potential energy. It is obvious that the temporal range of such a mechanism has a wide range of dependency. The propagation processes and the cause-and-effect relationship between various pre-seismic phenomena at different altitudes are still not very understandable. Therefore, it could be a possible reason for our observable parameters to have different pre-seismic time domains. Depending on the most sensitive and convincing parameter, one can have a sufficient idea of such a precursory time frame. For this, numerous numbers of EQs have to be studied through this process. In addition, some other significant parameters like ULF emissions, ozone concentration, atmospheric conductivity profile, etc. need to be studied additionally. This will improve the concept of the LAIC mechanism and give a better idea of solving this well-known problem of precursory phenomena of seismic hazards.

The importance and difficulties of the numerical model associated with the pre-seismic mechanism have already been mentioned. Not every physical model for LAIC has an identical hypothesis and produces similar results. As mentioned in the Introduction, Piersanti et al. [40] give highly convincing outcomes by using GPS-TEC and AGW variations. Though our outcomes for TEC and AGW are pre-seismic in nature, they differ from their findings. This is a similar problem of “difference in the time frame” and the fundamental driven force responsible for the change in TEC and AGW. Our work has a limitation of having only two IGS stations, and we are unable to generate a spatial distribution of such TEC variation. Of course, it will be a high priority to run our outcomes through their model for some other earthquake which will be done in the future. Secondly, as mentioned above, this work will also provide an opportunity to understand the internal mechanism of earthquake processes and how they get migrated to the different layers of the atmosphere producing different temporal variations for different seismogenic parameters. Thus, it will bring great motivation to understand the pathology of the inner earth, and thus this internal physical mechanism needs to be understood thoroughly. We will apply all such processes in the near future.