Analysis of Potential Precursory Pattern at Earth Surface and the Above Atmosphere and Ionosphere Preceding Two Mw ≥ 7 Earthquakes in Mexico in 2020–2021

Comprehensive analyses of seismic precursory characteristics may provide evidence for the possible mechanisms of lithosphere-atmosphere-ionosphere coupling (LAIC). The integrated monitoring data from space and ground measurements makes it possible to reveal the potential earthquake precursory pattern obeying LAIC. In this paper, we study seismic precursors from air temperature data of the Center for Environmental Prediction (NCEP), Global Navigation Satellite System (GNSS) based Total Electron Content (TEC) and electron density retrieved from Swarm satellite (Alpha) data in two Mw 𝐴𝐴 ≥ 7.0 earthquakes in Mexico in 2020–2021. We find seismic thermal anomalies can be identified from temperature variation in the lithosphere and atmosphere by improved tidal force fluctuant analysis (TFFA). The spatial-temporal evolution of seismic thermal anomalies is significantly different from that of non-seismic thermal anomalies. In addition, synergistic observations of the Swarm Alpha satellite and GNSS stations can provide more comprehensive information about perturbations in the ionosphere. The results show that anomalies originate along fault lineaments near the epicenter, leading to air pressure gradient and ionization process in the atmosphere, followed by ionospheric perturbations. Furthermore, composite anomalies of LAIC first appear around a month before earthquakes and tend to recur several times over earthquake preparation phases. Our study provides more evidence for understanding the characteristics of earthquake precursors. Key Points: • Spatial-temporal evolution of seismic thermal anomalies in lithosphere and atmosphere is investigated based on the variation of daily tidal force potential • The temporal rule of anomalous transmission supports the coupling of lithosphere-atmosphere-ionosphere • The earthquake precursory pattern obeying lithosphere-atmosphere-ionosphere coupling is explored

propose appropriate methods for detecting thermal anomalies associated with earthquakes (Chen et al., 2020;Tramutoli et al., 2005;Tronin et al., 2002). Recently, several researchers have applied tidal force fluctuant analysis (TFFA) to detect seismic thermal anomalies in the atmosphere Xu et al., 2021;Zhang et al., 2021), based on the indicative function of the tidal force in the seismic time-domain (Heaton, 1975). The methods can be further extended to the observation of seismic thermal anomalies in the lithosphere, but there are still some limitations to be resolved (Xu et al., 2021;Zhang et al., 2021).
The ionospheric perturbations during earthquakes preparation phases have been studied by different researchers. Shah et al. performed a detailed statistical analysis of Total Electron Content (TEC) anomalies before and after 1182 (Mw > 5.0) earthquakes during 1998-2019 based on Global Navigation Satellite System (GNSS) observations . The peak plasma frequency of ionosonde data was applied to three earthquakes in Pakistan and Nepal by Ahmed et al., to determine the duration of ionospheric anomalies (Ahmed et al., 2018). In addition, Parrot studied ion density data over the epicenters during nighttime from low altitude satellite DEMETER (Parrot, 2011), and proved the feasibility of detecting disturbances in the ionosphere through satellites. Similarly, Marchetti and Akhoondzadeh detected anomalies with parameters of the ionosphere from the Swarm satellites before a powerful earthquake with Mw = 8.2 and opened up the possibility of creating a monitoring system based on satellite data . Moreover, De Santis et al. reported the statistical significance of ionospheric anomalies observed by Swarm satellites as anticipating strong earthquakes on a global scale during 2014(De Santis, Marchetti, et al., 2019. The application of diverse observational data is significant for the identification of phenomena at different scales. The appearance of the LAIC may be related to some exchanges of mass and energy in the Earth system during earthquake preparation, which are often manifested by various and weak phenomena . Several works apply multiparametric analysis to explain the interaction between different phenom ena and further explore possible chains of physical phenomena before earthquakes Marchetti et al., 2020;Sasmal et al., 2021). There is no agreement on the evolutionary mechanism of LAIC, but the origination of LAIC is consistent across different theoretical reports, that is, fault activity Freund, 2011;Kuo et al., 2014;Pulinets & Ouzounov, 2011;Pulinets et al., 2018). The fault activity leads to the accumulation of energy from the squeezed rock under tectonic stress in seismogenic zones (Ma et al., 2008). Under these circumstances, convective thermal fluids are decomposed and rise through geological structures as heat flux, thermal water vapor, and gaseous emissions, eventually reaching the surface (Piscini et al., 2017;Pulinets & Ouzounov, 2011). The direct manifestation of fault activity during earthquakes preparation phases is the increase in the surface temperature (Ouzounov et al., 2006). However, the occurrence time of pre-seismic thermal anomalies is not consistent with the time of observed LAIC phenomenon in previous studies. Researchers have conducted several long-term statistical studies in different regions, and most pre-seismic thermal anomalies occurred more than a week before the mainshock days (Eleftheriou et al., 2016;Genzano et al., 2021;Zhang & Meng, 2019). In contrast, many LAIC phenomena occur more frequently with the approaching of earthquakes, and peak 5 days or less before earthquakes (Adil et al., 2021;Kumar et al., 2021;Mehdi et al., 2021;Shah, Aibar, et al., 2020;Yang & Hayakawa, 2020).
The predisposing factors of thermal anomalies are diverse during earthquakes preparation phases. Seismic activities always lead to synergistic variations of surface temperature rising and atmospheric warming, which start at the ground surface in the vicinity of different tectonic faults within the earthquake preparation area. The local influence has a characteristic evolution which is different from the temperature change caused by atmospheric advection. Seismic thermal anomalies can be identified more accurately according to the characteristic in spatial-temporal evolution. Hence, the synchronous criteria of the seismic thermal anomaly in the lithosphere and atmosphere deserves our attention. Furthermore, inadequate observation of ionospheric perturbations may limit the detection of LAIC. Space measurements are up to date and cover a wide area, while ground stations measurements have excellent time continuity and stability. The combination of observational data from different sources can provide a more explicit understanding of the possibilities involved in these phenomena. For more specific detection of ionospheric perturbations, the synchronized observation combining space and ground measurements should be considered.
In this paper, we explore the potential earthquake precursory pattern obeying LAIC linked to two Mw≥ 7.0 earthquakes in Mexico. To identify accurately seismic thermal anomalies in the lithosphere and atmosphere as successive air temperature variation, we apply TFFA and resolve its limitations. The combination of GNSS stations observations and Swarm satellite data is applied to the detection of disturbances in the ionosphere. Furthermore, we evaluate the path of anomalous transmission from the lithosphere to the ionosphere and the characteristic of the earthquake precursory pattern based on seismic thermal anomalies and ionospheric perturbations.

Seismic Data
As one of the most seismically active regions in the world, Mexico is located in a region with five regularly interacting plates: Cocos, Caribbean, Rivera, North America, and Pacific, see Figure 1a. The Cocos plate and Rivera Plate move northeastward, the Caribbean plate moves eastward, the North American plate moves southwestward, and the Eastern Pacific plate moves northwestward. Cocos plate subducts beneath the North American plate, and the subduction interface indicates a significant variety in strike dip angles along the central portion of the Middle American Trench (Perez-Oregon et al., 2021). Such interaction around the Mexican Pacific's southern coast is the main inducement of earthquakes in the country. The earthquakes studied here include an Mw = 7.4 earthquake occurred on 23 June 2020 and an Mw = 7.0 earthquake occurred on 8 September 2021, which were caused by reverse faulting and shallow thrust faulting near the plate boundary between Cocos and North American plates respectively, as shown in Figure 1b.
In this study, thermal anomalies and ionospheric perturbations related to two earthquakes (Mw ≥ 7.0) in Mexico are investigated during 3 tidal force periods. Specifically, the time interval considered for Mw = 7.4 earthquake is from 35 days before the earthquake to 7 days after the earthquake, and for Mw = 7.0 earthquake is from 38 days before the earthquake to 6 days after the earthquake. The information on earthquakes is obtained from the United States Geological Survey (https://earthquake.usgs.gov/earthquakes). At first, the method proposed by Genzano et al. (2021) is used to avoid earthquake clustering. This method is based on calculating the corresponding Dobrovolsky region (i.e., R = 10 0.43 ) of each earthquake (Dobrovolsky et al., 1979), as follows: • Dobrovolsky radius of each earthquake is calculated according to magnitude. • Earthquakes are sorted by increasing time of occurrence and decreasing magnitude. • Starting from the first event, the earthquakes that occurred within a month and fall into the corresponding Dobrovolsky radius are excluded. • This process is repeated until all the earthquakes have been evaluated.

Air Temperature
The air temperature images are acquired from National Center for Environmental Prediction (NCEP) reanalysis data, which is the product conducted by NCEP and the National Center for Atmospheric Research (NCAR) since 1991 (Kalnay et al., 1996). The NCEP reanalysis data are developed by assimilating various observations from Global Satellite Observation Data, Ocean Data, Global Ground Observation Data, Unmanned Aerial Vehicle Data, etc. The NCEP reanalysis data are available at the ground surface and 26 isobaric surfaces with 1° × 1° horizontal spatial resolution and 6-hr temporal resolution (Kistler et al., 2001). The ground surface and five isobaric surfaces are selected to detect the spatial and temporal anomalies of air temperature. Due to the various altitudes of epicenters, five isobaric surfaces with the shortest distance above the epicenter are extracted for each earthquake, as shown in Table 1.

Tidal Force and Thermal Anomalies
In this study, we explore the surface thermal anomalies and the atmospheric thermal anomalies. The atmosphere is subject to many short-term disturbances, which may cause fluctuations in air temperature. Therefore, it is necessary to adopt a method that can eliminate continuous spatiotemporal variation in the atmosphere. The occurrence of earthquakes originates from the mechanical processes that start with tectonic stress accumulation by subjecting the rocks inside the earth to reaching a critical breaking point. As an important stress parameter to trigger earthquakes, the tidal force induced by celestial bodies can cause tectonic stress in the rocks to achieve a critical state (Kilston & Knopoff, 1983). Also, the tidal force exerts an inducing effect on the formation of thermal anomalies during earthquake preparation by changing the state of seismotectonic stress . While affecting the tectonic stress build-up, the tidal force causes the squeezed rock to release energy, which in turn acts on the decomposition of convective thermal fluids. Tidal Force Fluctuant Analysis (TFFA) is a method to identify thermal anomalies based on tidal force periods. The procedure for calculating tidal force potential has been described in detail by Ma et al. and Zhang et al. (Ma et al., 2012;Zhang et al., 2018). TFFA divides data into different periods according to the variation of daily tidal force potential Xu et al., 2021;Zhang et al., 2021). However, the application of TFFA is limited by the need to determine the period from the relationship between earthquake occurrence time and tidal force potential (Xu et al., 2021). To use TFFA for earthquake precursor identification, we make the following simplification. Each top of tidal force is taken as the first day of the next period, and the day before the next top is taken as the last day of each period. For each position, the mean value of the whole period is set as the background. The difference between the temperature of single days and the temperature of the background is used to detect thermal anomalies. The specific criteria (Xu et al., 2021) for identifying thermal anomalies related to earthquakes are as follows: • Thermal anomaly emerges bottom-up. As height rises, the thermal anomaly tends to diminish. • The distribution of thermal anomaly at each vertical level coincides with the fault pattern. • Thermal anomaly spreads on the time scale, thereby thermal anomaly should last two consecutive days at least.

Total Electron Content (TEC)
In this study, TEC is retrieved from permanent stations of the International GNSS Service (IGS) (http://igs.org/ network). The 30 s temporal TEC from IGS permanent stations is obtained, which is then transformed into diurnal variation for subsequent analysis. TEC is described in TEC Units (TECU) with 1 TECU = 10 16 electron/m 2 . The computation of line-of-sight TEC (also called slant TEC or STEC) from raw GPS data is described by Dautermann et al. (2009) and Mannucci et al. (1999). The STEC is then transformed into vertical TEC (VTEC) as described by Mehdi  Furthermore, the confidence bounds based on the median ( ) and interquartile range (IQR ) are constructed to determine anomalies of TEC. The variation of TEC retrieved from the GNSS stations is extremely sensitive.
Despite the diurnal periodicity of TEC, there are often significant numerical differences in the segment with the same trend. Therefore, wider confidence bounds should be considered for TEC anomaly identification. The confidence bounds for a single day are calculated based on the data of 10 days before and 10 days after the observed day, as defined in Equation 1.
where UB is the upper bound, and is the lower bound. The difference between observed TEC and confidence bounds is calculated as deviations in TEC (dTEC) (Shah, Tariq, Ahmad, et al., 2019), as described in Equation 2.
Then, in order to further exclude the influence of subtle disturbances on anomaly identification, we extract the outliers of dTEC through standard deviation ( ). 2 is set as the upper bound, and −2 is set as the lower bound.

Electron Density From Swarm
The Swarm mission is composed of three identical quasi-polar orbiting satellites (Alpha, Bravo, and Charlie, also named Swarm-A, Swarm-B, and Swarm-C, respectively) (Friis-Christensen et al., 2006). Specifically, Swarm-A and Swarm-C fly almost parallel at lower orbit (around 460 km altitude), and Swarm-B flies at higher orbit (around 510 km altitude). The electron density data used in this study are measured by the Langmuir probe (LP) of Swarm-A with a sampling rate of 2 Hz. For each earthquake, the sample points of the track of Swarm-A within the Dobrovolsky region are extracted from the data. Then, the quality flags are used to judge whether the electron density data of sampling points are useable or not, and all the useable data are sorted out. The differences between electron density measured by Swarm-A and the corresponding values modeled by IRI2016 (Bilitza et al., 2017), defined as ∆ E, are calculated. The median values of ∆ E separately for day and night tracks are calculated for each day. Finally, the time series of median values for day and night are constructed. To detect the anomalies of ∆ E, the median ( ) and the interquartile range ( ) are applied. The time series consists of the median values of ∆ E, indicating that the influence of transient disturbances has been ruled out. Thus, the confidence bounds should be narrowed to detect abnormal electron density measurements, as described in Equation 3 where UB is the upper bound, and LB is the lower bound.

Results
We investigate various anomalies in the lithosphere, atmosphere, and ionosphere during three periods of each earthquake, and the evidence of the potential earthquake precursory pattern has been revealed. The studied periods of two earthquakes in Mexico are divided according to variation of daily tidal force potential, as shown in Figure 2. The temperature variation during the studied periods is shown in Figures S1 and S2 of Supporting Information S1. According to specific criteria (Xu et al., 2021), the thermal anomalies related to the two earthquakes are extracted. Figures 3 and 4 show the seismic thermal anomalies that occurred during the studied periods for Mw = 7.0 earthquake (8 September 2021) and Mw = 7.4 earthquake (23 June 2020), respectively. For the Mw = 7.0 earthquake, there are four seismic thermal anomalies with various duration (from Aug 1 to 3, from Aug 11 to 13, from Aug 16 to 21, and from Aug 24 to 27). For the Mw = 7.4 earthquake, there are three seismic thermal anomalies with various duration (from May 19 to 28, from June 2 to 10, and from June 17 to 23). It is observed that most seismic thermal anomalies follow the progressive change rule of "initial rise-intensification-attenuation," which  is corresponding to the thermal infrared radiation procession of rocks that broke under stress loading (Wu et al., 2006). This phenomenon is also noted by Ma et al. (2018). The absence of the rule in other seismic thermal anomalies may be caused by two reasons: (a) the seismic thermal anomalies exist for a short time, so the process cannot be reflected; (b) the process of the seismic thermal anomalies is covered by other stronger anomalies in the same period. Furthermore, non-seismic thermal anomalies are observed in two forms (see Figures S1 and S2 in Supporting Information S1): (a) In period C of the Mw = 7.0 earthquake and periods A and B of the Mw = 7.4 earthquake, thermal anomalies appear in the atmosphere and decay up to down. This form of thermal anomaly, caused by a strong meteorological phenomenon, exists for a long time and tends to weaken gradually; (b) In periods A and C of the Mw = 7.0 earthquake, we can see the non-seismic thermal anomalies that initiate on the surface. This form of thermal anomaly exists for a short time, does not distribute along the active faults, and changes irregularly.
We have implemented the identification of seismic thermal anomalies in the lithosphere and atmosphere. To verify LAIC, abnormal variations of TEC and electron density are detected. Previous studies show that the enhancement of TEC related to earthquake is more prominent over the earthquake breeding zones Shah & Jin, 2015;Xie et al., 2021). In addition, Xie et al. further reported that TEC anomalies were positively correlated with electron density anomalies (Xie et al., 2021). Therefore, in this study, the positive anomalies of TEC and electron density are extracted and used for subsequent analysis. As seen in Figure 5, TEC anomalies appear at the beginning of the study period for the Mw = 7.0 earthquake, and the geomagnetic condition remains quiet that day. The peak of dTEC appears 12 days before the earthquake. Unfortunately, it is dubious the association of the disturbance with earthquake due to highly disturbed geomagnetic conditions (Dst less than −70 and Kp greater than 4). Similarly, a disturbance in the ionosphere that appears at early phases is observed in the Mw = 7.4 earthquake ( Figure 6). Then, the anomalies repeat two times at both GUAT station and MANA station under quiet geomagnetic condition. It is noted that the anomalies observed at MANA station, which is further away from the epicenter, are weaker and have a slight delay compared to those observed at GUAT station. A similar association between the distance from the epicenter and ionospheric perturbation is noted by Kumar et al. (2021). Figures 7 and 8 illustrate the Swarm-A day-time and night-time electron density residual variations during studied periods of the two earthquakes. A striking anomaly is detected at night 37 days before the Mw = 7.0 earthquake ( Figure 7). However, the appearance of this anomaly may be related to disturbed geomagnetic conditions (Dst greater than 35 and Kp greater than 4). Then, the ∆ E increases significantly on three single days (28, 19, and 6 days before the earthquake) and two consecutive days (from 22 to 21 days before the earthquake), when geomagnetic activities remain quiet. The first electron density anomaly occurs at night 27 days before the Mw = 7.4 earthquake (Figure 8), and then the ∆ E shows an upward trend on three single days (27, 21, and 16 days before the earthquake) and two consecutive days (from 4 to 3 days before the earthquake).
To provide more evidence to the analysis of the earthquake precursor pattern obeying LAIC, seismic thermal anomalies and positive anomalies of TEC and electron density are compared on time scales, as shown in Figure 9. It should be noted that only TEC anomalies occurring at both GNSS stations at similar times are recorded for the Mw = 7.4 earthquake, and the day with high geomagnetic activity is marked in magenta. The statistics in Figure 9 presents that most anomalies of TEC and electron density fall within the range of seismic thermal anomalies. Meanwhile, one TEC anomaly and one electron density anomaly fall outside the range of seismic thermal anomalies for Mw = 7.0 earthquake (Figure 9a). Data from only one GNSS station are available for detecting the anomalies of the Mw = 7.0 earthquake, and this may lead to the result that local deviations in the ionosphere caused by meteorological activities are observed (Ahmed et al., 2018;Rishbeth, 2006;Rishbeth et al., 2009). This local deviation is also reflected at the stations used to observe the Mw = 7.4 earthquake. Specifically, an anomaly is observed at MANA station 7 days before the earthquake but not at GUAT station at a similar time. The electron density anomaly outside the range of seismic thermal anomalies occurred 6 days before the earthquake (the third day of period C). There are non-seismic thermal anomalies around that day, which could cause the seismic thermal anomalies to be covered.

Discussion
We have investigated seismic and non-seismic thermal anomalies in the lithosphere and atmosphere during the preparation phases of two Mw≥ 7.0 earthquakes in Mexico. Compared with non-seismic thermal anomalies, seismic thermal anomalies have significant discrepancies in spatial-temporal evolution. It should be noted that the seismic thermal anomalies observed recur over multiple tidal force periods. In previous studies, most seismic thermal anomalies occurred only in a single tidal force period before the earthquake Xu et al., 2021). The process of seismic thermal anomalies follows the rule of radiation of rock fracture under stress loading . During the preparation of an earthquake, the process by which stress builds up leading to rock fracture, resulting in energy release, is reduplicative, thereby seismic thermal anomalies usually occur several times before earthquakes (Fu et al., 2020;Panda et al., 2007;Tronin et al., 2002). The original TFFA considers the first day of a tidal force period as the background, making the detection susceptible to meteorological activity, whereas TFFA applied in this study reduces the influence of meteorological factors by averaging the tidal force period. Thus, improved TFFA can observe the occurrence of seismic thermal anomalies over multiple tidal force periods. Furthermore, the limitation of extracting anomalies only after earthquakes have occurred is also resolved in TFFA.
In addition, another work in this study includes detecting ionospheric perturbations associated with earthquakes. GNSS stations and Swarm-A satellite data are used to observe variations of TEC and electron density, respectively. Due to the impact of geomagnetic storms on the ionosphere and Earth magnetic field (Afraimovich & Astafyeva, 2008), diurnal variation of geomagnetic indices before and after every earthquake is investigated. According to limit values of geomagnetic indices, the TEC anomaly 12 days before Mw = 7.0 earthquake and the electron density anomaly 37 days before Mw = 7.0 earthquake may be related to geomagnetic activity, while other abnormal changes of TEC and electron density are not caused by geomagnetic activity. The TEC anomalies occurred for the first time about a month before the two earthquakes. The disturbances being detected such early are consistent with past studies (  the disturbances may have occurred much earlier . The last occurrences of TEC anomalies precede that of electron density anomalies, and the occurrence time of TEC anomalies and electron density anomalies is inconsistent. This is related to the observation method of the two parameters. GNSS station observes single point data with continuous time, while Swarm-A satellite observes partial time data with a continuous track. Further analysis is based on the integration of thermal anomalies and ionospheric perturbations. The comparison on the time scale is presented in Figure 9. The anomalies of TEC and electron density outside the range of seismic thermal anomalies may be caused by local deviation and seismic thermal anomaly being covered up, respectively (Figure 9a), as described in Section 3. It should be noted that most disturbances in the ionosphere are observed during seismic thermal anomalies. This temporal relationship provokes thoughts about the path of anomalous transmission and the earthquake precursory pattern. The generation of ionospheric perturbations lag behind that of seismic thermal anomalies, and ionospheric perturbations disappear before seismic thermal anomalies. This indicates that there is a mechanism affecting the ionosphere in the formation of seismic thermal anomalies. The anomalies initiate in the lithosphere and their distribution does not depend on the epicenter territory location, which is consistent with the lithosphere-atmosphere-ionosphere coupling (LAIC) (Pulinets & Ouzounov, 2011). The distribution of seismic thermal anomalies along fault lineaments corresponds to the tectonic stress accumulation and loss during the preparation phases (Sykes et al., 1999). Seismic thermal anomalies in the atmosphere evolve from surface due to thermal radiation of rock fracture under the stress loading. Xu et al. (2021) reported the coupling of "temperature rise," "relative humidity decrease" and "air mass uplift" in the atmosphere before earthquakes based on TFFA (Xu et al., 2021). Such coupling is also a characteristic of the ionization process in the atmosphere (Pulinets et al., 2018). This indicates that the ionospheric perturbations before earthquakes could be caused by the air pressure gradient (Genzano et al., 2007), which in turn is induced by seismic thermal anomalies in the atmosphere. On the other hand, we observe reduplicative seismic thermal anomalies accompanied by ionospheric perturbations starting around a month before earthquakes. The occurrence time of the first precursors obeying LAIC in this study is consistent with the first appearance of surface thermal anomalies in the past statistics (Eleftheriou et al., 2016;Genzano et al., 2021;Zhang & Meng, 2019). Moreover, the repeatability of earthquake precursors initiates from thermal anomalies induced by fault activity. This shows that repeated observations of earthquake precursors can be performed over a longer time scale. Certainly, multi-source data, appropriate methods, and clear anomaly recognition criteria are necessary. Plenty of impending anomalies are observed in previous case studies (Cicerone et al., 2009;Daneshvar & Freund, 2017;Genzano et al., 2007;Mehdi et al., 2021;Shah, Aibar, et al., 2020), whereas the signals are sometimes strong, more often fleeting (Freund, 2011). The earthquake precursory pattern with the characteristics including LAIC, repeatability, and early occurrence is significant for earthquake prediction.

Conclusions
In this paper, we have investigated the appropriate methods to detect anomalous patterns related to earthquakes covering the lithosphere, atmosphere, and ionosphere during the seismic preparation period. The results show that modified TFFA reduces the influence of meteorological factors so that seismic thermal anomalies can be observed over multiple tidal force periods. In addition, the simplified period division method makes it possible to detect earthquake precursors using TFFA. In ionospheric observation, different from using a single data source,  perturbations can be identified more accurately by combining the temporal continuity of GNSS station data with the spatial continuity of satellite data. Under the action of thermal radiation from rock fracture, seismic thermal anomalies initiate along the active faults lines and affect various atmospheric parameters, and then influence the ionosphere. Thus, seismic thermal anomalies are strongly correlated with ionospheric perturbations, indicating the following anomalous transmission path: thermal anomalies originate along fault lineaments near epicenter and are then transmitted to atmosphere, leading to generation of air pressure gradient and ionization process, and followed by ionospheric perturbations. Further analysis demonstrates that there is an interesting earthquake precursory pattern. The occurrences of ionospheric perturbations observed by GNSS stations and satellites are well covered by the duration of seismic thermal anomalies. Such composite anomalies obeying LAIC first appear around a month before earthquakes and tend to recur several times over earthquake preparation phases. This pattern provides more feasibility for earthquake precursor identification and earthquake prediction. This study puts forward an important step moving from the observation of LAIC phenomenon to the search for the earthquake precursory pattern obeying LAIC. The data and methods applied in our experiment are the solution proposed here, but other tactics may resolve better and more efficiently the same goal. Further studies still need to be undertaken on other large earthquakes.

Data Availability Statement
The air temperature data after tidal force fluctuant analysis used for thermal anomaly detection, the TEC data and electron density data used for detection of ionospheric perturbations in this study are available at Dryad via https://doi.org/10.5061/dryad.37pvmcvmv with CC0 1.0 Universal (CC0 1.0) Public Domain Dedication license (Xu, 2022).