Seismogenic Anomalies in Atmospheric Gravity Waves as Observed from SABER/TIMED Satellite during Large Earthquakes

Department of Ionospheric Sciences, Indian Centre for Space Physics, 43, Chalantika, Garia Station Road, Kolkata 700084, India Department of Electrical and Electronics Engineering, Ancient Olive Grove Campus, University of West Attica, 12241 Egaleo, Greece Department of Space Science and Engineering, National Central University, Taoyuan 32001, Taiwan Hayakawa Institute of Seismo Electromagnetics, Co. Ltd, UEC Alliance Center, Chofu Tokyo 182-0026, Japan


Introduction
According to several research, it is well accepted that there are significant disturbances detected in the regular atmospheric processes due to seismic hazards [1][2][3]. A hypothetical mechanism called Lithosphere Atmosphere Ionosphere Coupling (LAIC) has been proposed to investigate and validate this strange physical mechanism that exists underneath the earth's surface and propagates up to ionospheric heights through a number of geophysical and geochemical processes. LAIC is proposed to be functional through three major channels: (a) the chemical, (b) the acoustic, and (c) the electromagnetic [3]. There are a wide variety of processes involving electromagnetic disturbances, starting from disruption in the very low frequency (VLF)/low frequency radio signals and ultralow frequency (ULF)/extreme low frequency (ELF) [4][5][6][7][8][9][10][11][12], changes in plasma density in higher ionospheric altitudes [13][14][15], irregularities in total electron content (TEC) [16][17][18][19][20][21][22] etc. The primary acting agent of the acoustic channel is atmospheric gravity wave (AGW) excitation, which may occur as an outcome of atmospheric oscillation in stratospheric heights above the epicenter of a similar earthquake, swaying upward, and disrupting upper atmospheric altitudes. Gravity waves (GWs) are generated as the gravitational force or buoyancy attempts to reestablish harmony in a liquid medium or at the interface of two mediums, according to fluid dynamics [23,24]. GWs are the method for transmitting energy to the stratosphere and mesosphere from the lower atmosphere. Two of the principal factors in GW generation are frontal systems or the wind stream moving through mountains in the lower atmosphere. The momentum transfer function is as follows: first, waves spread with an almost continuous mean velocity across the atmosphere. The wave amplitude increases as you rise higher in height because of the low air density. Nonlinear force breaks the waves and transfers the momentum to the medium [25]. This energy exchange is responsible for driving some of the atmosphere's largest dynamical highlights. This GW affects stratospheric climate, temperature, and pressing factors in the lower stratosphere and is energized by convection systems, fly currents, and fronts that vary in size between a few hundred meters to a few kilometers, with a duration that falls between the Brunt Vaisala period and the inertial period [25]. Profound convective clouds can also create AGWs because of the collaboration between unsteady convective movements and the encompassing stable surroundings [25]. These waves change in scale from that of individual convective updrafts to bigger scopes characterized by coordinated convective structures [26]. Convectively created GWs can influence the middle atmosphere's momentum, can produce turbulence and mix, and can connect with the general climate to advance or smother new convection [27,28]. The factors such as sudden stratospheric warming and meteorological events can also generate the AGW. AGW oscillations, which arise around the earthquake epicenter, cause fluctuations in temperature, pressure, ground motion, and other factors to disrupt the atmosphere. AGWs have the ability to propagate upward, causing ionosphere disruptions [7,20]. AGW is a dispersion branch in atmospheric waves with a period of around 5 minutes to 10 hours and a wavelength ranging from~10 meters to 1000 kilometers. Due to the viscous dissipation of the short-wave components, its wavelength increases with altitude, reaching hundreds of meters in the ionospheric D layer (between~60 and~90 kilometers above sea level) and~10 kilometers in the F2 layer (above~400 kilometers). These are the fastest oscillations in the atmosphere, capable of generating disturbances that can reach the ionosphere [29,30].

Related Works
Several studies on the preseismic activity of AGW have already been reported over the last few decades. This section briefly describes those works. Before an earthquake, the hypothesis of atmospheric wave excitation due to seismic activity is reported by Garmash et al. [31], Linkov et al. [32], and Shalimov [33]. Numerous studies using both satellite and ground-based observations were conducted to examine the preseismic AGW activity. Between 1985 and 1989, AGW energies were investigated using the middle and upper atmosphere (MU) radar. Additionally, the observational results demonstrate that the jet stream exhibits an annual variance [34]. Later, Tsuda et al. [35] summarized the characteristics and variations of GW with height, season, and latitude from a combination of the abovementioned MU radar, two medium frequency (MF) radars, and a lidar. Tsuda et al. [36] investigated the global distribution of potential energy across mid-latitude using a global positioning system (GPS) temperature profile and observed that potential energy is higher during winter seasons. Dhaka et al. [37,38] have used the Mesosphere-Stratosphere-Troposphere (MST) radar at Gadanki, India, to report AGW activity and its relationship with convections during Indian southwest rainfall. Miyaki et al. [39] were the first to have reported the AGW signature from subionospheric low-frequency signals and its role in the LAIC mechanism in Japan. Korepanov et al. [40] concluded that AGW can be an important parameter in the seismoionospheric study by using data on surface atmospheric pressure and magnetic field. To examine AGW-related potential energy ðE p Þ, Zhang et al. [41] analyzed the temperature data from the Sounding of the Atmosphere Using Broadband Emission Radiometry (SABER) instrument installed on the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) satellite. Nakamura et al. [42] have performed a comparative investigation and attempted to track down the relating seismogenic effect for a few earthquakes. For the 2004 Niigata-Chuetsu earthquake (M = 6:8), wavelet analysis of those parameters reveals modifications over a period of 10 to 100 minutes, which is consistent with AGW. Using nighttime VLF fluctuation and the SABER temperature profile, previously AGW activities were investigated prior to the Imphal earthquake. Around a week before this earthquake, a large increase in E p linked to the AGW was seen near the epicenter [43]. Piersanti et al. [20] observed enhanced AGW activity on the day of the 2018 Bayan earthquake computed from ERA5 data and also found anomalous activity in GPS-TEC on the same day. Carbone et al. [29] used a mathematical model lithosphere-atmosphere coupling for seismic events to compute vertical atmospheric temperature profile and compared the obtained results with observation. Sasmal et al. [22] found abnormal AGW activity 6 days before the 2020 Samos earthquake using the SABER temperature profile and 11 days prior to the earthquake using wavelet analysis of GNSS-TEC. Previously, Ouzounov et al. [44] used ground-and space-based observation to investigate the preseismic anomaly. The majority of the research are based on the LAIC mechanism. The data from outgoing longwave radiation (OLR), GPS-TEC, and foF2 are employed in the research. In OLR data, a large amount of radiation was emitted on March 7, an anomalous increase in TEC was seen on March 8, and abnormalities in f0F2 were recorded between March 03 and 11, 2011. For the 2011 Tohoku earthquake, Ohta et al. [45] used the Chubu University ULF/ELF network to investigate ULF/ELF electromagnetic radiation. On March 6, 2011, before the earthquake, ULF/ELF atmospheric radiation was detected. Two days before this earthquake, Sasmal et al. [46] noticed a significant shift in the sunrise  15,2016), and Chiapas earthquake (September 08, 2017) in various channels of LAIC mechanism. So, it is interesting to see the effect of these abovementioned earthquakes in the acoustic channel. This paper presents a detailed analysis of AGW activity in the stratosphere using SABER temperature profile for the specific earthquakes. The detailed analysis of this work is presented in Section 2. The findings of this research are presented in Section 3. Finally, the concluding remarks are presented in Section 4.

Materials and Methods
According to United States Geological Survey (USGS) (https://earthquake.usgs.gov/), the earthquakes are classified into seven different classes according to their magnitude. Based on this classification, the earthquake magnitudes (M) between 7.0 and 7.9 are classified to the "major earthquake" class, while M ≥ 8 earthquakes are classified to the "great earthquake" class. Therefore, the studied 2011 Tohoku earthquake (M = 9), 2012 Indian Ocean earthquakes (M = 8:6 and 8.3), and 2017 Chiapas earthquake are great earthquakes, while Nepal earthquakes (M = 7:8 and 7.3) and Kumamoto earthquake (M = 7) are major earthquakes. The abovementioned earthquakes were chosen to study the AGW activity. The general information of the earthquakes as collected from USGS are presented in Table 1. Each earthquake's seismic preparation zone was calculated using Dobrovolsky's radius equation. The equation provides the radius of a circle representing the preparation zone of the earthquake as: R = 10 0:43M , where M is the magnitude of the earthquake. The epicenter of the earthquake is situated at the center of the circle [56]. The preparation zone of each one of the studied earthquakes is presented in Table 1 and shown in Figure 1.
The To avoid possible contamination of the atmosphere due to solar geomagnetic storms, the geomagnetic situation is studied for ±30 days from the day of earthquake. D st and K P are the main geomagnetic indices during such period 3 Journal of Sensors for individual earthquake is shown in Figure 2. The geomagnetic data (Kp and, D st ) are taken from OMNIWEB NASA archive (https://omniweb.gsfc.nasa.gov/). Figures 2(a) and 2(g) show that three mild geomagnetic storms were detected during the Tohoku earthquake on February 14, March 1, and March 11, 2011. The first geomagnetic storm occurred 24 days before the earthquake with a minimum D st value of -40 nT and a maximum Kp value of 5.3. The earthquake occurred 10 days after the second geomagnetic storm. For the specific geomagnetic storm, the D st has a minimum value of -88 nT, and the maximum value of Kp is 5.3 for this storm. On the earthquake day, there was a third geomagnetic storm. D st lowest value was -83 nT, and the highest value of Kp was 5.3. Figure 2(b) shows that the D st reached a minimum value of -60 nT on April 5, 2012, which is six days before the earthquake and a second minimum value of -60 nT on April 13, 2012, which is two days after the earthquake. As seen in Figure 2(h), the maximum value of Kp is 4.0 and 4.3, which is less than 5.0, indicating that both storms are minor. Therefore, the preearthquake period was quiet for the Indian Ocean earthquake. Again, the minimum value of D st for the storm that happened 14 days after the earthquake was -120 nT. For this storm, the maximum value of K P is 5.7. Figures 2(c) and 2(i) and Figures 2(d) and 2(j) reveal that the seven days prior to both earthquakes were geomagnetically quiet for the Nepal earthquake.     [67,68]. There are several measures involved in extracting AGW from a temperature profile. The temperature (T) profile for a specific geographic area with a finite latitude/longitude range is first retrieved from the SABER archive data (http://saber.gats-inc.com/). The logarithm of each individual temperature profile is then computed. Fitted values are derived from the logarithmic temperature using a third-degree polynomial. The residual values are then derived by subtracting the fitted temperature profile from the initial profile. A 4 km boxcar filter is added to the residual values to remove other waves of wavelength shorter than 4 km. The final profile is computed by combining the filtered residual values with the fitted value. The least square fit (LSF) of these profiles is the antilogarithm. The regular zonal mean temperature and other zonal wave components are calculated using these LSFs, ranging from 1 to 5. The number of all wave components from 0 to 5 is the background temperature (T 0 ). By subtracting the background temperature from the initial profiles, the perturbation temperature (T p ) is obtained. By putting the values into equation (1), the potential energy (E p ) associated with the AGW is derived.  Journal of Sensors where g denotes gravity's acceleration, and N denotes the Brunt-Väisälä frequency described by Here, z is the altitude, and c p is the specific heat at constant pressure.
A nine-dimensional matrix is generated after the computation of E p . Latitude, longitude, date or day of the year in UT, altitude, original SABER temperature profile, reconstructed fitted temperature profile, perturbation        11 Journal of Sensors used to determine background changes. Three years preceding the earthquake year were used for this calculation. At various altitudes, the recorded E p values differ. The average of E p corresponding to the three years was used to determine the E p background threshold value. If the E p value of any day in the ± 15 days period from the day of the earthquake exceeds the threshold value, then this day is considered to be an anomalous day possibly associated with the specific earthquake. It is worth noting that any anomalous day observed by this method must have been geomagnetically and meteorologically quiet. All the metalogical conditions around the earthquake period were gathered from Japan Meteorological Agency (JMA)     13 Journal of Sensors day is identified on April 5, 2016. The maximum and threshold value of E p for the specific earthquake is 7.2 J/kg and 3.06 J/kg, respectively, at 42 km altitude. The highest E p value for the Chiapas earthquake is about 9.3 J/kg at 42 km on August 25, 2017, with the threshold value 5.6 J/kg. All the E p variations for the identified anomalous days are presented in Figure 3.
The time altitude variation of E p associated with AGW is computed from the nine-dimensional matrix for the abovementioned earthquakes, and the three-dimensional profile of such E p variation is presented in Figures 4-9. The X -axis shows the dates (in UT) around the earthquakes, and Y-axis represents the altitude profiles in km the color bars show the E p in J/kg. For the Tohoku earthquake, Figure 9 shows a significant increment in E p that occurred around 43 km from March 5 to 7, 2011, that is approximately 4 to 6 days prior to the earthquake. A significant increase of E p occurred after the earthquake day (March 20-23, 2011). Another significant increment is observed from February 12 to 14 around 39 to 41 km. This increment of E p may be associated with the earthquake (M > 5) which occurred on February 16 in the same region as confirmed by USGS. Figure 5 shows a significant enhancement of E p at around 34 to 39 km from April 5 to 9, 2012, which occurs 2 to 6 days   (May 12, 2015). For the Kumamoto earthquake, Figure 10 shows that the E p is anomalous increased around 40 to 42 km from April 4 to 6, 2016 about one week prior to the earthquake. For the Chiapas earthquake, there is continuous enhancement occurring before the earthquake. Figure 8 shows that enhancement of E p is significant around 45 to 48 km during August 17 to 23, 2017. The next enhanced activity is observed around 40 to 42 km from August 23 to 25, 2017. Another enhancement occurs around 41 to 45 km from August 29, 2017, to September 2, 2017. The same effect of equatorial KWs contaminates those seismogenic AGW excitation [69] for the Chiapas earthquake, as it does for the Indian Ocean earthquake (Figure 6).

Spatial variation of AGWs.
After detection of the altitude of maximum the E p , the main focus is on the spatial variation of E p over the epicentral region for all the earthquakes. The spatial distributions of E p around the Tohoku  For Kumamoto earthquake, it is observed that E p is maximum at the epicentral region on April 6, April 11, and April 15, 2016, at an altitude of 43 km which reflects from the spatial distribution of Figure 18 during April 1 to 15, 2016. For the Chiapas earthquake, Figure 19 shows that the E p reaches its maximum value on August 19, 2017, at 41 km altitude and after that, it started to diminish at the epicentral region; again, it is observed on August 25, 2017, around from west direction and scattered around east direction on the next day. The same kind of result is found at 42 km on the same days which are shown in Figure 20. From Figure 21, it is evident that the enhancements in E p is maximum during August 17 to 18, 21 to 25, 28, and 30 to 31, 2017 at 46 km near the epicentral region.

Discussions
In this manuscript, the study is mainly focused on preseismic activity of AGW before four great earthquakes and three major earthquakes. It is evident from Figures 4 and 8 Figure 22. Convective systems usually exist around cyclones causing AGWs. In a similar way, for the Kumamoto earthquake, the enhancement of E p on the earthquake day is not generated due to the geomagnetic storm that occurred on the earthquake day. It is mainly formed due to a few frontal systems accompanied cyclones passed the studied region (20°-30°N, 120°-140°E) during April 14 to 19, 2016 shown in the weather maps ( Figure 23).
In spatial variation of AGW, the huge patches seen in Figures 11-21 are mostly the result of synoptic and planetary-scale (hundreds to thousands of kilometers) meteorological systems. These systems can cause disturbances in

18
Journal of Sensors the atmosphere and either move fast eastward or westward. As a result, they generate wide-scale temperature fluctuations, resulting in huge patches of high E p region, as seen in these figures.

Conclusions
As reported by various researchers in recent years, the possibility of AGWs playing a significant role in detecting preseismic disturbances is indeed a well-proven phenomenon [7,8]. AGW is one of the regulating agents in the acoustics channel of the LAIC mechanism, which originates from pressure or temperature convection caused by seismogenic sources. This manuscript elaborately discusses such phenomena based on a space-based observation. An indirect approach is used to detect this AGW operation from the principle of convective currents throughout the lower stratosphere since there is no clear way to classify the excitation of AGW in the stratosphere before earthquakes. SABER/ TIMED temperature profile is used to calculate the E p associated with such AGW in this manuscript. We focused on four great and three major earthquakes occurred at various periods and in various geographical regions around the world. The spatiotemporal patterns of potential energy are investigated for these seven earthquakes that occurred in Tohoku, the Indian Ocean, Nepal, Kumamoto, and Chiapas.
For the Tohoku earthquake (March 11, 2011), AGW excitation in a form of potential energy enhancement is observed at 43 km altitude and on 4 to 6 days prior to the earthquake day. For the Indian Ocean earthquake (April 11, 2012), similar AGW activity became maximum on 2 to 6 days prior to the days of the earthquake at 34 to 37 km altitude. For the Nepal earthquake, the preseismic irregularities started a bit earlier. The first AGW activity has been detected on 13 to 14 days before the day of the first Nepal earthquake (April 25, 2015) at 35 to 36 km altitude, and another enhancement has been identified on 1 to 3 days prior to the day of the same earthquake at an altitude 40 to 42 km. For the second Nepal earthquake (May 12, 2015), an enhancement of potential energy associated with AGW has been observed on 3 to 4 days prior the earthquake day at an altitude of 34 to 36 km. For the Kumamoto earthquake (April 15, 2016), at a height of 43 km, AGW activity has been detected, on 8 to 9 days prior to the earthquake day. At an altitude of 42 to 46 km, the anomalous increase in E p correlated with AGW was observed one week prior to the Chiapas earthquake (September 8, 2017). Some postearthquake AGW activity has been observed for the Tohoku and the Kumamoto earthquakes. These anomalies due to typhoons and cyclones over the Japan region are reported during the same time frame and responsible to generate such wavelike phenomena. The equatorial KWs also play a vital role in AGW excitation [69] Significant enhancement in E p has

19
Journal of Sensors been observed after the main shock for the Indian Ocean earthquake as it happened near the equatorial region. For the Chiapas earthquake, a similar contribution from equatorial KWs possibly contaminated the result.
According to JMA, no cyclones or convective systems were formed in the 15 days prior to the 2011 Tohoku earthquake. From March 21 to 23, many cyclonic systems emerged at 45°-60°N in the postearthquake period. In this study, for Tohoku earthquake, the maximum enhancement of E p is observed before the earthquake around March 05 to 07, 2011. As no cyclones occurred during that time frame, the increased AGW activity may have been caused by the earthquake. The postenhancement of E p between March 21 and 23 was caused by the formation of cyclonic systems in that region. There were no cyclonic systems formed in the preearthquake timeframe, according to Bureau of Meteorology reports of March and April 2012. From April 16 to 25, 2012, a postearthquake low pressure area (classified as 19 U) was developed in Australia, with the lowest pressure of 1006 hPa. In our observation, the highest enhancement of E p is observed before the earthquake around April 5 to 9, 2012 and since the period is meteorologically and geomagnetically quiet, the generation of AGW is possibly connected to the earthquake. The postenhanced activity of AGW may be connected to that low pressure formation around Australia or may be due to contamination of KWs and aftershocks of the earthquake. According to weather report of April and May from IMD, no cyclones were formed between the chosen pre-and postearthquake period around that region for both Nepal earthquakes. Since the period was meteorologically as well as geomagnetically quiet, the enhanced activity of AGW around April 22-24 and May 09, 2015, may be connected to earthquake. According to the JMA weather report for April 2016, no cyclonic systems were formed before the earthquake; however, cyclonic activity was reported in the selected region from March 14 to 19,2016. From the observation, the enhanced AGW activity is observed around April 4 to 6, 2016. The activity may be linked to the earthquake because the time is free of all kinds of meteorological phenomena and geomagnetically quiet. The other enhancements around April 14 to 19 may be connected to the cyclonic system around the region. According to NWS and NOAA reports, the preearthquake period was meteorologically quiet and free from any kind of cyclone or convective systems but a major hurricane Irma passed around 20°N to 25°N and -70°W to 85°W within spatial area chosen for Chiapas earthquake from September 07 to 09, 2017. Another major hurricane Katia also passed through the chosen area around 22°N to 25°N and 94°W to 98°W from September 05 to 09, 2017. Though the periods August 17-19 and August 31-September 2 were meteorologically quiet with solar activities, Sasmal et al. [22] and Biswas et al. [43] have demonstrated that the geomagnetic storm has a contaminating impact in the electromagnetic channel of LAIC but not in the acoustic channel, implying that AGW activity may be unaffected. So, the observed activity is possibly connected to the earthquake. The period from August 24 to 28 was absolutely solar quiet; so, the enhanced activity of AGW is possibly connected to earthquake.
The observed outcomes are compared with the previously detected AGWs for the same earthquake. The study of Yang et al. [51] used ERA5 reanalysis data for the investigation of AGW activity for the Kumamoto earthquake and presented a similar increase in AGW activity 4 to 6 days prior to the earthquake. In SABER observation, similar kind of results both in the temporal and spatial measurements is identified. In another work by Yang et al. [49] for 2011, the Tohoku earthquake is also verified by the observation through SABER. Wave-like structure in the range of AGW has been detected previously in a different way by Chakraborty et al. [9] where the computation of AGW has been achieved by taking the wavelet of nighttime VLF signal during Nepal earthquake. A significant enhancement of AGW 3 days before the second Nepal earthquake (May 12, 2015) is observed in Chakraborty et al. [9]. Very recently, Biswas et al. [43] also compared the AGW enhancement both observed from night-time VLF fluctuation as well as SABER observation for a separate earthquake in Imphal, India, on January 3, 2016 (UT) where the identical spatiotemporal effect has been unidentified. Thus, the observation of this study through a new methodology is consistent with the previous works and, for most of the cases, the AGW enhancement is preseismic in nature. As the physical mechanism behind the LAIC is still under investigation, this study needs much attention. In contrast to earlier works, this works dealt with a large number of earthquakes for investigating such AGW phenomena. As, in most cases, there are no such huge space weather phenomena that happened before the earthquakes, and these AGW activities are mostly due to the seismogenic sources. To testify these preseismic AGW variations as a key ingredient of the acoustics channel of LAIC, other parameters from thermal (surface latent heat flux, relative humidity, etc.) and electromagnetic channel (VLF radio anomalies, total electron content (TEC) anomalies, energetic particle precipitation, etc.) are being thoroughly analyzed separately and will be reported elsewhere as a multiparametric study of LAIC phenomena.

Data Availability
The data of atmospheric temperature is taken from the SABER data archive (http://saber.gats-inc.com/), the geomagnetic data is taken from OMNIWEB NASA Archive, and the earthquake regarding information is taken from USGS.

Disclosure
This paper is focused on SABER satellite observations of the precursory effect for massive earthquakes in stratospheric atmospheric gravity wave (AGW).