Atmospheric Anomaly Analysis Related to Ms > 6.0 Earthquakes in China during 2020–2021

: The attention towards links of atmospheric parameter variation and earthquakes has increased exponentially by utilizing new methods and more accurate observations. Persistent research makes it possible to gain insight into the precursor mechanism of earthquakes. In this paper, we studied the universality of detecting atmospheric anomalies associated with earthquakes based on tidal force ﬂuctuation in China for earthquakes of Ms > 6.0, and explored the inﬂuence of tidal force on tectonic stress. The data of air temperature, geopotential height, ozone mixing ratio, and relative humidity from the National Center for Environmental Prediction (NCEP) were analyzed to reveal the spatiotemporal variation of atmospheric anomalies at multiple isobaric surfaces. Furthermore, the coupling of atmospheric parameters was investigated. The results showed that continuous solicitation exerted by tidal forces could change the strength of tectonic stress that causes earthquakes. The evolution pattern of air temperature, geopotential height, and relative humidity could be supported by atmospheric thermal vertical diffusion, while the anomalies of ozone mixing ratio was not evident. This veriﬁed the feasibility of detecting multi-parameter atmospheric anomalies associated with earthquakes based on tidal force ﬂuctuation. Our results provide more evidence for understanding the atmospheric precursor characteristics of earthquakes. atmospheric anomalies, also showing that measurement may play an important role for the analysis and prediction of Further studies are recommended to reveal the atmospheric anomalies related to with multiple atmospheric parameters and remote sensing


Introduction
There are numerous anomaly reports associated with earthquakes, not only on earth surface but multiple disturbances are also reported in atmosphere during seismic preparation period [1][2][3][4][5][6][7][8][9][10]. In particular, thermal anomaly has been widely studied as the impending earthquake anomaly because of its correlation with tectonic activity and high temporal resolution of thermal infrared satellite data [11]. Various theoretical reports were conducted to propose the methods for identifying thermal anomalies related to earthquakes. In the early research of seismic thermal anomaly detection, different methods for interpreting large areas were applied to infrared and land surface temperature measurements [2,3,12,13]. The Robust Satellite Technique (RST) data analysis approach was developed based on the multi-temporal analysis of a historical dataset to detect seismic thermal anomalies [14,15]. Wavelet transform was also applied to extract seismic thermal anomalies by researchers [16,17]. However, the above methods depend on historical data of a long time series, thereby the measurement of small thermal variation related to earthquake is considered unreliably [18]. Furthermore, the images acquired by thermal infrared sensors onboard satellites are affected by high cloud cover [19][20][21]. Although detection of thermal anomalies is often disrupted by historical data and cloud, these obstacles can be avoided through proper methodology and background data. The tidal forces induced by celestial bodies can cause tectonic stress in the rocks to achieve critical state that triggers earthquakes [22]. The changes in tidal forces and thermal anomaly phenomenon related to seismic activity reflect the mutation of tectonic activities [18,23]. A related research shows that tidal forces, as the method of pre-calculating Earth deformation phenomena, serves

Tectonic Environment of Earthquakes
In the study, a total of six earthquakes (Ms > 6.0) in China during 2020-2021 are investigated with atmospheric measurements during 5 periods of tidal force potential. The information of the earthquakes is obtained from the China Seismographic Network and United States Geological Survey. At first, the method proposed by Genzano et al. [35] was used for earthquake declustering. The geographical location of epicenters is shown in Figure 1. Meanwhile, the detailed information of earthquakes is listed in Table 1.

Meteorological Data
The atmospheric measurement data are provided by the National Center for Environmental Prediction (NCEP)as a product of atmospheric data conducted by NCEP and the National Center for Atmospheric Research (NCAR) since 1991 [36]. The NCEP data are developed after assimilating meteorological observations from the global satellite observation data, global ground observation data, ocean data, radiosonde data, and unmanned aerial vehicle data, etc. The NCEP data, available at 28 pressure levels with 1 • × 1 • horizontal spatial resolution and 6-hour temporal resolution, offer the most comprehensive, global reanalysis of meteorological data [37].
In this study, spatial and temporal anomalies of air temperature, geopotential height, ozone mixing ratio, and relative humidity are investigated over the earthquake preparation zone to explore atmospheric anomalies induced by earthquake. In order to further investigate variability in the atmosphere caused by earthquakes, multiple layers of atmospheric parameters are selected. Due to the various altitude of different epicenters, five isobaric surfaces with the shortest distance above the epicenter were selected for each earthquake, as shown in Table 2.

Meteorological Data
The atmospheric measurement data are provided by the National Center for Environmental Prediction (NCEP)as a product of atmospheric data conducted by NCEP and the National Center for Atmospheric Research (NCAR) since 1991 [36]. The NCEP data are developed after assimilating meteorological observations from the global satellite observation data, global ground observation data, ocean data, radiosonde data, and unmanned aerial vehicle data, etc. The NCEP data, available at 28 pressure levels with 1°×1° horizontal spatial resolution and 6-hour temporal resolution, offer the most comprehensive, global reanalysis of meteorological data [37].
In this study, spatial and temporal anomalies of air temperature, geopotential height, ozone mixing ratio, and relative humidity are investigated over the earthquake preparation zone to explore atmospheric anomalies induced by earthquake. In order to further investigate variability in the atmosphere caused by earthquakes, multiple layers of atmospheric parameters are selected. Due to the various altitude of different epicenters, five isobaric surfaces with the shortest distance above the epicenter were selected for each earthquake, as shown in Table 2.

Tidal Force and the TFFA Method
The tidal force induced by the sun and the moon have macroscopic effects on the earth, and it is one of the important external causes of ground stress [38]. Tidal force causes periodic stress changes in the earth's interior. The accuracy description of calculating the tidal force potential can be found in [23,39]. In the study, the tidal force potential of epicenters was considered at a time scale of 24 h, and the calculation time of daily tidal force potential was consistent with earthquake occurrence time, as shown in Table 1.
The TFFA method proposed by Ma et al. [39] was applied to extract thermal anomalies related to earthquakes, and it was described in detail by Zhang et al [25]. TFFA first divides the data into different periods according to the variation of daily tidal force potential. The first day of each period is set as the background, and thermal anomalies associated with Remote Sens. 2021, 13, 4052 4 of 23 earthquakes are identified by comparing other data with the background. However, the method can be further modified in terms of period division. In the modified method, the top or trough adjacent to the day of earthquake is considered as the first day of each period. Genzano et al. analyzed thermal anomalies of earthquakes in Japan for 11 consecutive years. They found that few thermal anomalies occurred near the day of earthquakes occurrence, and most earthquakes were apart from thermal anomalies at least a week [35]. Therefore, the trough or top near earthquake is more suitable as the background. The modified period division of TFFA is described as:

•
If the earthquake occurred at the top of tidal force, then each top is taken as the first day of each period. Similarly, if the earthquake occurred at the trough, then each trough is taken as the first day of each period.

•
If the earthquake did not occur at the trough or the top, time distance is used as the standard. We choose the one closer to the time of earthquake. If the trough and top are equidistant from the time of earthquake, the preceding trough or top is taken as the first day of each period.
After period division, the variations of atmospheric parameters can be calculated based on the background. ∆T, ∆G, ∆H, and ∆O are defined as Equation (1), and the first day of each period is considered as background.
where T d is the temperature of a day except background, T background is the temperature of the background, G d is the geopotential height of a day except background, G background is the geopotential height of the background, H d is the relative humidity of a day except background, H background is the relative humidity of the background, O d is the ozone mixing ratio of a day except background, and O background is the ozone mixing ratio of the background.

Specific Criteria for Selecting Atmospheric Thermal Anomalies Related to Earthquakes
Atmospheric temperature is affected by many factors, and it has many forms of evolution. In particular, Ma et al. [26] indicated that the process of atmospheric thermal anomalies related to earthquakes is similar to thermal infrared (TIR) variation of rocks that broke under stress loading. The identification of atmospheric thermal anomalies requires the definition of a threshold. The evolution of atmospheric thermal anomalies is a continuous process, in fact, thermal anomalies appear at the bottom and dissipate in the upper layer. Atmospheric thermal anomalies have a steep change in temperature from high to low on the spatio-temporal scale [26]. High threshold could make the evolution pattern ignored. The problem can be avoided by obtaining a set of atmospheric thermal anomalies with a lower threshold and then identifying the thermal anomalies associated with earthquakes according to the evolution pattern. Therefore, the threshold ϕ is set as Equation (2).
where ∆T is the mean of ∆T for the entire study area within the total study duration, and σ ∆T is the standard deviation of ∆T for the entire study area within the total study duration. In order to identify atmospheric thermal anomalies related to earthquakes, we established the specific criteria including spatial distribution, temporal persistence, and relative intensity:

•
The threshold ϕ is used to detect thermal anomalies; • Thermal anomalies appear at the bottom and lift towards the top. As height increases, the thermal anomalies diminish; Remote Sens. 2021, 13, 4052 5 of 23 • The distribution of thermal anomaly at each vertical level coincide with the fault pattern; • Thermal anomalies should be extended on time scale, so that thermal anomalies should last two consecutive days at least.
Atmospheric thermal anomalies of each earthquake were identified according to above criteria, and then ∆G, ∆H, and ∆O of corresponding days were also listed to analyze atmospheric evolution.

Change in Tidal Force and Determination of Period
The tidal force potential of epicenters was calculated at a time scale of 24 h, and the calculation of tidal force potential was performed at earthquake occurrence time every day. The daily tidal force potential at the hypocenter of studied earthquakes was calculated, as shown in Figure 2. For each earthquake, the time series was divided into five periods, and each period was labeled as A, B, C, D, and E, respectively. The Ms 7.5 earthquake in Maduo and the Ms 6.1 earthquake in Biru occurred at the inflection points of period D, and other earthquakes occurred around the day when the tidal force potential was at peak. In addition, the specific dates of each period for earthquakes are listed in Table 3.

Atmosphere Evolution during Atmospheric Thermal Anomalies
The specific criteria for selecting atmospheric thermal anomalies related to earthquakes was applied to all periods of each earthquake, and atmospheric thermal anomalies were identified in four earthquakes, as shown in Table 4. In earthquakes of Maduo and Yutian, no atmospheric thermal anomaly was detected. In Yangbi, the last day of atmospheric thermal anomaly occurred 30 days before the earthquake. In Biru, the last day of atmospheric thermal anomaly occurred 19 days before the earthquake. In Nima, the last day of atmospheric thermal anomaly occurred 33 days before the earthquake. In Jiashi, the last day of atmospheric thermal anomaly occurred 21 days before the earthquake.
In order to identify atmospheric thermal anomalies related to earthquakes, we established the specific criteria including spatial distribution, temporal persistence, and relative intensity:  The threshold φ is used to detect thermal anomalies;  Thermal anomalies appear at the bottom and lift towards the top. As height increases, the thermal anomalies diminish;  The distribution of thermal anomaly at each vertical level coincide with the fault pattern;  Thermal anomalies should be extended on time scale, so that thermal anomalies should last two consecutive days at least. Atmospheric thermal anomalies of each earthquake were identified according to above criteria, and then ∆ , ∆ , and ∆ of corresponding days were also listed to analyze atmospheric evolution.

Change in Tidal Force and Determination of Period
The tidal force potential of epicenters was calculated at a time scale of 24 h, and the calculation of tidal force potential was performed at earthquake occurrence time every day. The daily tidal force potential at the hypocenter of studied earthquakes was calculated, as shown in Figure 2. For each earthquake, the time series was divided into five periods, and each period was labeled as A, B, C, D, and E, respectively. The Ms 7.5 earthquake in Maduo and the Ms 6.1 earthquake in Biru occurred at the inflection points of period D, and other earthquakes occurred around the day when the tidal force potential was at peak. In addition, the specific dates of each period for earthquakes are listed in Table 3.

Atmosphere Evolution During Atmospheric Thermal Anomalies
The specific criteria for selecting atmospheric thermal anomalies related to earthquakes was applied to all periods of each earthquake, and atmospheric thermal anomalies were identified in four earthquakes, as shown in Table 4. In earthquakes of Maduo and Yutian, no atmospheric thermal anomaly was detected. In Yangbi, the last day of atmospheric thermal anomaly occurred 30 days before the earthquake. In Biru, the last day of atmospheric thermal anomaly occurred 19 days before the earthquake. In Nima, the last day of atmospheric thermal anomaly occurred 33 days before the earthquake. In Jiashi, the last day of atmospheric thermal anomaly occurred 21 days before the earthquake.   In order to provide more evidence to atmospheric anomalies, other processes within the atmosphere over each earthquake were monitored by spatial and temporal analysis of geopotential height, ozone mixing ratio and relative humidity during atmospheric thermal anomalies. The most obvious thermal anomalies and corresponding atmospheric parameters are shown in Figures 3-6, and the composite atmospheric evolution was pre-sented in Figures A1-A4 (dynamic evolution of atmospheric parameters is in Appendix A). Figures 3-6 correspond to earthquakes in Yangbi, Biru, Nima, and Jiashi, respectively.

Discussion
Earthquakes occur due to the mechanical process that begins with tectonic stress by subjecting the rocks inside the earth to accumulating energy, culminating in a destructive rupture. Tidal forces act as an external factor, which supplements tectonic stress, exerting During atmospheric thermal anomalies, their distribution was consistent with fault pattern for each earthquake. The atmospheric thermal anomalies in Yangbi lasted throughout period A. It is obvious that thermal anomalies appeared at the bottom and lift towards the top with a fading trend. For geopotential height, the trend was not evident on 7 April. Since 8 April, the geopotential height on the top isobaric surface had been abnormal in the position similar to thermal anomalies. The trend of geopotential height anomalies was opposite to that of thermal anomalies, which gradually disappears from the top to the bottom. This phenomenon could be observed throughout the duration of thermal anomalies. In the period, the variation of ozone mixing ratio had no observable anomalies. In contrast, the negative anomalies of relative humidity over the epicenter were consistent with the thermal anomalies on time scale and spatial distribution. Similar geopotential height and relative humidity anomalies were also evident during thermal anomalies in Biru, and the variance of ozone mixing ratio showed no correlation with the earthquake. In Biru, the opposite trend of geopotential height compared with thermal anomalies was also evident, and geopotential height anomalies peaked on 26 February, one day after the peak of thermal anomalies. In addition, the trend of geopotential height in Nima and Jiashi accorded with that in Yangbi and Biru, while the trend of relative humidity in only Jiashi accorded with that in Yangbi and Biru. The relative humidity anomalies in Nima was not evident in upper layers, and the relative humidity anomalies had a weak similarity with thermal anomalies in spatial distribution. The ozone mixing ratio anomaly was also not observed in Nima and Jiashi.
In summary, the distribution of geopotential height and relative humidity was observed to be similar to that one of thermal anomalies, but the distribution of ozone mixing ratio did not show consistency with that one of thermal anomalies. The positive anomalies of geopotential height had an obvious opposite trend to thermal anomalies, and the trend of negative relative humidity anomalies was consistent with that one of thermal anomalies. In addition, a delay of geopotential height anomalies relative to thermal anomalies was observed.

Discussion
Earthquakes occur due to the mechanical process that begins with tectonic stress by subjecting the rocks inside the earth to accumulating energy, culminating in a destructive rupture. Tidal forces act as an external factor, which supplements tectonic stress, exerting an obvious inducing effect on the active fault whose ground stress was at the critical state [23]. Among the studied earthquakes, four of them occurred around peaks of tidal forces, and two of them occurred at the inflection points of tidal forces. This reflects that although tidal forces have an inducing effect on earthquakes, they cannot trigger earthquake directly. The distribution of atmospheric thermal anomalies in Yangbi, Biru, Nima, and Jiashi obeyed the rule of thermal radiation of rock broken under the stress loading, indicating that atmospheric thermal anomalies are possibly produced by the preparation phase of earthquakes. The atmospheric thermal anomalies of earthquakes in Yangbi and Nima occurred on the first day after the peaks of tidal forces, while the atmospheric thermal anomalies of earthquakes in Biru and Jiashi occurred near the troughs of tidal forces. The continuous solicitation of tidal forces and accumulation of tectonic stress could lead to rock micro-fracture and earthquakes [26]. Therefore, rock micro-fracture and earthquakes are more likely to occur around peaks of tidal forces, but not always.
In addition, the variation of other atmospheric parameters (i.e., geopotential height, ozone mixing ratio, and relative humidity) were observed at isobaric surfaces during atmospheric thermal anomalies. It can be seen that all atmospheric thermal anomalies related to earthquakes were accompanied by consistent geopotential height variation, that is, the positive anomalies of geopotential height appeared at the highest isobaric surface and decayed downward. Heating of the lower isobaric surface is favorable for air mass raising, which can lead to the increasing of geopotential height at upper isobaric surface [40,41]. The atmospheric thermal vertical diffusion associated with the earthquakes involves air being warmed by land, uplifted by heat flux, cooled and dissipated in the sky [26]. The geopotential height variation pattern observed during atmospheric thermal anomalies accords with atmospheric thermal vertical diffusion and air mass migration. In contrast, the relationship between ozone mixing ratio and thermal anomalies was not observed. The content of tropospheric ozone is much lower compared with that of stratospheric ozone, which makes it difficult to distinguish the abnormal variation of tropospheric ozone. The occurrence of tropospheric ozone variance depends on many environmental factors (e.g. dust storms [42], volcanic eruptions [43], etc.), thereby the previous monitoring of tropospheric ozone is based on comparative analysis of long time-series data [44][45][46]. In particular, the research by Piscini et al. shows that the total column ozone was less statistically significant through analyzing long time-series data [44], and our study indicates that the anomalies of tropospheric ozone in short time-series are also not evident. Therefore, tropospheric ozone may be not an effective parameter that precedes the earthquake occurrence. The negative anomalies of relative humidity observed in Yangbi, Biru, and Jiashi were consistent with thermal anomalies on time scale and spatial distribution, whereas the negative anomalies of relative humidity observed in Nima were observed on the similar time scale as thermal anomalies. Heating of land promoted air movement, air mixing, and an increase in air temperature. This causes the ions to attract water vapor molecules, leading to a decrease in the content of free available water vapor [4,47]. As a result of the phenomenon, negative anomalies of relative humidity are generated.
It should be noted that atmospheric parameters are influenced by many external factors. Strong meteorological phenomena, geological activities, and disasters can lead to significant variation of atmospheric parameters [42,43,48,49]. The air temperature at the top isobaric surface in Nima rises again in part of the time. Such short-term changes may be related to some meteorological activities, and the weak anomalies of relative humidity in Nima may also be caused by this. Furthermore, the absence of observed thermal anomalies in Maduo and Yutian may be also associated with meteorological activities. During the observation of Maduo (3 April−17 June), the weather is undergoing seasonal variation. Due to latitude and altitude, seasonal weather changes greatly increased air temperature in Maduo. In June 2020, a severe heat wave was reported in Yutian. Therefore, the obvious seasonal climate change and short-term intense climate change could affect the detection of atmospheric thermal anomalies. Furthermore, the earthquakes in the study were all shallow (depth ≤ 10 km). The anomalies associated with the earthquakes are supposed to be influenced by depth [13]; therefore, further studies are necessary to reveal the influences that depth of earthquakes changes on TFFA.
Among the six studied earthquakes, the atmospheric thermal anomalies were observed at four locations. This suggest that it is feasible to identify atmospheric thermal anomalies based on tidal force fluctuation. Furthermore, the coupling relationship of thermal anomalies, geopotential height and relative humidity was also observed, indicating that TFFA can be further applied to the identification of atmospheric anomalies associated with earthquakes. However, the detection of atmospheric anomalies is affected by many environmental factors, and the anomalies may be completely covered. The determination of periods in TFFA is also limited by the need to know when an earthquake occurred. This limits the use of the method in forwarding direction to identify earthquake precursors, but the method is significant in the post-analysis of earthquake precursors.

Conclusions
In this paper, earthquake-induced atmospheric anomalies have been detected from air temperature, geopotential height, ozone mixing ratio, and relative humidity for six earthquakes based on tidal force fluctuation. Research has shown that it is feasible to identify atmospheric thermal anomalies using TFFA, which is not only suitable for specific earthquakes. Tidal forces can induce rock micro-fracture and earthquakes, but they cannot trigger rock micro-fracture and earthquakes directly. The continuous solicitation exerted by tidal forces can reduce the strength of tectonic stress that causes earthquakes. Correspondingly, the atmospheric thermal anomalies detected in four seismogenic zones only occurred in one of five periods. This implies that tidal forces cause atmosphere temperature changes by altering stress state rather than directly triggering thermal anomalies. Further analysis has demonstrated the coupling relationship of thermal anomalies, geopotential height, and relative humidity, and no anomaly of tropospheric ozone was observed in this study. Atmospheric thermal anomalies obeyed the rule of thermal radiation of rock broken under the stress loading at several locations, providing more evidence for the tectonic stress accumulation and loss in the seismogenic zone. In addition, the positive anomaly variation pattern of geopotential height and the negative anomaly variation pattern of relative humidity could be supported by atmospheric thermal vertical diffusion, and tropospheric ozone may be not an effective parameter that precedes the earthquake occurrence. Coupling phenomena of atmospheric parameters were detected in several seismogenic zones, which could help to identify earthquake atmospheric anomalies, also showing that atmospheric measurement may play an important role for the analysis and prediction of earthquakes. Further studies are recommended to reveal the atmospheric anomalies related to earthquakes with multiple atmospheric parameters and remote sensing data.  Figure A1. Cont.