Temperature Variations in Multiple Air Layers before the Mw 6.2 2014 Ludian Earthquake, Yunnan, China

: On 3 August 2014, an Mw 6.2 earthquake occurred in Ludian, Yunnan Province, China (27.245 ◦ N 103.427 ◦ E). This damaging earthquake caused approximately 400 fatalities, 1800 injuries, and the destruction of at least 12,000 houses. Using air temperature data of the National Center for Environmental Prediction (NCEP) and the tidal force ﬂuctuant analysis (TFFA) method, we derive the temperature variations in multiple air layers between before and after the Ludian earthquake. In the spatial range of 30 ◦ × 30 ◦ (12 ◦ –42 ◦ N, 88 ◦ –118 ◦ E) of China, a thermal anomaly appeared only on or near the epicenter before earthquake, and air was heated from the land, then uplifted by a heat ﬂux, and then cooled and dissipated upon rising. With the approaching earthquake, the duration and range of the thermal anomaly during each tidal cycle was found to increase, and the amplitude of the thermal anomaly varied with the tidal force potential: air temperature was found to rise during the negative phase of the tidal force potential, to reach peak at its trough, and to attenuate when the tidal force potential was rising again. A signiﬁcance test supports the hypothesis that the thermal anomalies are physically related to Ludian earthquakes rather than being coincidences. Based on these results, we argue that the change of air temperature could reﬂect the stress changes modulated under the tidal force. Moreover, unlike the thermal infrared remote sensing data, the air temperature data provided by NCEP are not affected by clouds, so it has a clear advantage for monitoring the pre-earthquake temperature variation in cloudy areas.


Introduction
In the past thirty years, there have been numerous reports about surface temperature changes prior to earthquakes [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Using various thermal infrared (TIR) remote sensing data, researchers have proposed several methods for extracting thermal anomalies before earthquakes. For instance, Ouzounov and Freund [6] proposed the LST (land surface temperature) differencing method based on two years of data; Tramutoli et al. (2007) presented the robust satellite techniques (RST); Zhang et al. [15] used the wavelet transform method to extract the thermal infrared anomalies, and so on. Concerning TIR imagery data, images acquired by different thermal infrared sensors equipping different satellites have been widely used for the extraction of TIR anomalies before earthquakes, such as the Moderate-Resolution Imaging Spectroradiometer (MODIS) loaded on Terra and Aqua, the MeteoSat Second Generation-Spinning Enhanced Visible and Infrared Imager (MSG-SEVIRI), Advanced Very High Resolution Radiometer (AVHRR) loaded on the National Oceanic and Atmospheric Administration (NOAA) satellite, Chinese Feng Yun (FY) meteorological satellite, and so on [2,[16][17][18][19][20][21][22][23][24][25]. However, these methods are based on the statistical processing of multiyear remote sensing data, and they usually could not reliably measure small but significant thermal changes occurring prior to single earthquakes [26]. Blackett et al. [27] also found that the reported thermal infrared anomalies before the Gujarat (India) earthquake of 2001 were greatly influenced by cloud cover or other data gaps. Moreover, the TIR remote sensing is badly affected by clouds, and it is impossible to obtain the thermal information emitted by the ground when the interested area is cloudy. Zhang and Meng [28] found that the Sichuan Basin is often covered by clouds, and the seismic activity of this area is frequent, so that identifying thermal anomalies using thermal infrared data in such an area is limited.
On the other hand, the change in tidal force can trigger shallow thrust fault earthquakes [29], and the high correlation between earthquakes and tidal force is statistically significant and physically reasonable [30]. The tidal force will cause tectonic stress in rocks to reach the critical failure point [31][32][33], as the seismotectonic stress is the controlling factor for earthquake nucleation [26]. Some researchers have considered the variation of the tidal force to derive the time reference for choosing the background level above which to extract thermal anomalies, and then found that the extracted TIR anomalies and the tidal force were closely related to the seismic activity [26,[34][35][36][37]. Ma et al. [26] used multiple layers of air temperature to study the temperature change before, during, and after the Jiujiang earthquake, and found that seismic activity always leads to atmosphere warming before the earthquake, with heat originating from the ground, followed by atmosphere uplift diffusion, and finally dissipation in the sky.
In this paper, we first introduce in Section 2 the basic information about the Ludian earthquake. In Section 3, the variation of tidal force and the tidal force fluctuant analysis (TFFA) method are presented. The results of thermal anomalies before the Ludian earthquake and the long-term significance test for determining the relation between earthquakes and thermal anomalies are shown in Section 4, while the discussion and conclusion are presented in Sections 5 and 6.

Tectonic Environment of Ludian Earthquake
As shown in Figure 1, an Mw 6.2 earthquake occurred on 3 August 2014, 8:30:13 (UTC) in Yunnan Province, China, as a result of shallow strike-slip faulting within the crust of the Eurasia plate, with its epicenter located at 27.245 • N 103.427 • E. In the following 20 days, there were 18 aftershocks near Ludian and Yongshan with magnitudes ranging from 4.0 to 5.0. Tectonics here are broadly controlled by the convergence of the India and Eurasia plates, which has driven the uplift of the Himalayas to the west of this earthquake, and caused the formation of numerous intraplate continental transform structures in the surrounding region, including the left-lateral Xiaojiang fault system. Local deformation rates in the region of this earthquake indicate as much as 10 mm/yr of left-lateral shear, and at the latitude of the 3 August 2014 event, India moves northwestward with respect to Eurasia at a rate of approximately 51 mm/yr (https://earthquake.usgs.gov/earthquakes (accessed on 1 January 2021)).

Change of Tidal Force and the TFFA Method
The differential gravitational forces provided by the Sun and the Moon lead to deformations on the surface of the Earth, along its radial axis, by up to 40 cm [31], and the tidal force has a significant influence on the stress within the Earth. As shown in Figure 2, we calculated the hourly tidal force potential at the hypocenter of the Ludian event, and the space derivative of the potential is the tidal stress. More information about how to calculate the tidal force potential can be found in References [30,34]. Beeler and Lockner [38] studied the correlation between the failure times and oscillating stress based on laboratory experiments in which faults are loaded to quasiperiodic failure by the combined action of a constant stressing rate, intended to simulate tectonic loading, and a small sinusoidal stress, analogous to the tidal force. They found that when the period of the oscillating stress is less than t n , the correlation is not consistent with the threshold failure models, and the physical interpretation of t n is the duration of failure nucleation. Beeler and Lockner [38] indicated that "behavior at the higher frequencies is consistent with a second-order dependence of the fault strength on sliding rate, which determines the duration of nucleation and damps the response to stress change at frequencies greater than 1/t n ", and t n exceeds the daily tidal period. Therefore, rather than the hourly tidal force potential, we considered the tidal force potential at a time scale of 24 h to study the Ludian earthquake. Therefore, we selected the tidal force potential at 8:30:13 (UTC) every day; 8:30:13 (UTC) is the time when the Ludian earthquake happened. Ma et al. [34] proposed using the TFFA method to extract the thermal infrared (TIR) or thermal anomalies related to earthquakes. The first step of TFFA is to divide the data into different cycles according to the variation of daily tidal force potential. In TFFA, if the aiming earthquake occurred at the trough of sinusoidal tidal force, then each top of the sinusoidal variations is taken as the first day of each cycle. Similarly, if the earthquake occurred at the top, then the trough will be taken as the first day of each cycle. However, if the earthquake did not occur at the trough or the top, then the nearest preceding trough or top to the date of the earthquake is the candidate for the first day of each cycle. If the date of the earthquake is closer to the time of the trough rather than the top, we choose the trough, otherwise we choose the top. In this way, five fortnight cycles were based on the sinusoidal daily tidal force potential, as shown in Figure 2, each of them labeled as A, B, C, D, and E. As the red curve of Figure 2 shows, the Mw 6.2 earthquake occurred at the inflection point of cycle D. The five fortnight cycles we consider cover 16 June-1 July, 2-15 July, 16-30 July, 31 July-13 August, and 14-29 August. TFFA considers the first day of every cycle as the background reference, and ∆T values are the air temperature data or TIR data of other days subtracted by their corresponding background; then, the variation of ∆T is obtained. Finally, we need to set a threshold Φ to identify the TIR or thermal anomalies, if ∆T ≥ Φ, then it is regarded as an anomaly. The threshold is determined by the mean ∆T and standard deviation σ ∆T of ∆T for the entire study area within the total study duration, and Φ = ∆T + 3 * σ ∆T .

Multiple Layers of Air Temperature Evolution
The air temperature data are provided by NCEP, and it is a reanalysis product based on the Global Satellite Observation Data, Global Ground Observation Data, Integrated Ocean Data, Unmanned Aerial Vehicle (UAV) Observation Data, etc. The NCEP data are the most comprehensive, global reanalysis of meteorological data with 1 • × 1 • horizontal resolution, and this dataset provides temperature data for 26 isobaric surfaces [39]. Using the TFFA method, we can get the 3D variation of ∆T. The multiple layers can help us to distinguish thermal anomalies caused by earthquakes from the ones influenced by a change of weather. If the air temperature increases first in the upper layer and is much larger than in the lower layer, it is obvious that this rise is not caused by land heating and earthquakes, but by the local weather [26]. In this study, we selected five layers of air temperature to show the occurrence, uplift, dissipation, and disappearance of thermal anomalies along the vertical direction. The five layers are located, from lowest to highest, 2 m above the ground, and then at altitudes corresponding to 975, 925, 850, and 750 mb pressure. The reason why isobaric layers started at 975 mb is that the altitude of epicenter is approximately 350 m, which corresponds to the isobaric surface with 976 mb, and the layers above 2 m above the ground should be with lower pressure, so we choose 975 mb as the starting isobaric layer. Moreover, five layers of air temperature ending at 750 mb are sufficient for presenting the vertical occurrence-dissipation-disappearance of the anomalies, and in this case there will be no more anomalies in the layers with pressure lower than 750 mb.
First, as shown in the Figure 3, we do a simple statistical analysis to determine the threshold Φ to identify the anomalies. The negative skewness and positive kurtosis indicate that the distribution of ∆T is not Gaussian, and the Jarque-Bera test also rejects the hypothesis that the distribution of ∆T is Gaussian. In this case, it is shown in Figure 3 that Φ = ∆T + 3 * σ ∆T ≈ 11 K. Furthermore, as shown in Table 1, the proportion of ∆T values that are greater than or equal to 11 K is only about 0.57%. From a statistical view, in this case, these ∆T values that are larger than or equal to Φ (11 K) are obviously highly anomalous signals.

Description
Proportion of Qualified ∆T in Total ∆ As shown in Figure 4, there are thermal anomalies during cycles A, B, and C (as defined on Figure 2). All these anomalies appear in the lower layers, ascend, and then disappear in the high layers, so that they are identified as related to fault activity. We estimated the daily total precipitation near the epicenter (26 • -28 • N, 102 • -104 • E), which is shown in Figure 5 and at the bottom of Figure 6. Moreover, to quantitatively analyze the relation between thermal anomalies and other factors, we calculated the area of anomalies with ∆T ≥ 11 K in different layers, as shown in Figure 6. As shown in Figure 4 and at the top of Figure 6, the maximum area and duration (A: from 24 to 26 June, 3 days; B: from 4 to 13 July, 10 days; C: from 19 to 30 July, 12 days) of each cycle were gradually increasing (C > B > A) when the earthquake was approaching. In each cycle before the earthquake, air temperature began to rise when the tidal force was quickly decreasing, reaching a peak at the troughs of the tidal force, attenuated, and disappeared when the tidal force was recovering. In periods A and B, the thermal anomalies disappeared after the trough of tidal force; however, the thermal anomalies in period C did not; these had attenuated compared to the last day, and there was another round of rising temperature until the end of the period. Based on the trend, we speculate that there should be some stronger thermal anomalies on 31 July, 1, 2, and 3 August, because seismotectonic stress is severe and approaching the critical breaking point. However, there are no thermal anomalies in period D. As shown in Figure 5, there was heavy rainfall around the epicenter on 1 August, which caused severe landslides (https://earthquake.usgs.gov/earthquakes (accessed on 1 January 2021)). The bottom of Figure 6 also indicates that rainfall around the epicenter is not frequent in these three months, but the strongest rainfall in this period occurred on 1 August, and the amount of precipitation was much larger than that of other days. The heavy regional rainfall may thus have concealed the potentially strong thermal anomalies in period D. There were thermal anomalies in period E as well. However, as shown in Figure 4E, their spatial distribution was scattered, and they were not persistent, so they were not identified as related to earthquakes.
We cannot directly get the seismotectonic stress before the earthquake, so we choose the surface strain data provided by the China Earthquake Administration (CEA), which are obtained by the model YRY-4 four-component borehole strain gauge, to reflect the variation of stress before and after the earthquake. These data are not open access. The strain data in this case were recorded by the Zhaotong station (27.4 • N, 103.32 • E), which is the closest station to the hypocenter; the sensor lies 45 m underground and far from the depth of the hypocenter (12 km). As shown in the upper plot of Figure 6, the negative strain indicates that the observation point is under compression and contraction during this period; moreover, the negative enhancement of strain indicates that the compression grew more intense from 16 June to 29 August. Chen et al. [40] conducted an experiment to study the relation between stress-strain and variation of temperature, and the results showed that the surface strain is positively correlated to the strain strength. However, as shown in Figure 4 and at the top of Figure 6, there were no thermal anomalies after the earthquake when the strain was still negative and increasing. One possible hypothesis that explains this contradiction is that the deformation of surface ground (45 m) cannot reflect the variation of stress of deep layers near the hypocenters. That is, the elastic deformation of the surface ground is not the source of the thermal anomalies in this case; if it were, then there would be more thermal anomalies in periods D and E after the earthquakes.
Moreover, the variation of stress or strain around the hypocenter, where there might be a source creating the thermal anomalies, is different from the strain at the observation point (45 m). The mechanism of thermal anomalies created by earthquakes in the natural environment needs further study.

Significance Test for Thermal Anomalies
We introduce a significance test to prove whether or not the observed thermal anomalies are related to the Ludian earthquake. The variations of fortnight tidal force potential for different hypocenters and times differ much from each other. As shown in Figure 7, Figure 7 indicates that the backgrounds indicated by the TFFA method are spatially and temporally varying. This shows that the TFFA will be ineffective for areas far from the hypocenters and for different times when performing long-term statistical studies. Therefore, we selected the area around the epicenter (27 • -28 • N 103 • -104 • E) as the area of interest, and the average of ∆T for these four pixels from December 2007 to December 2014 was based on the TFFA method. We then set a relative intensity rule and a temporal persistence rule for extracting thermal anomalies: • Relative intensity rule-similar to the thermal anomalies of the Ludian earthquake, we set 11 K as the threshold; • Temporal persistence rule-unlike some thermal anomalies caused by short-term and instantaneous change of weather, the thermal anomalies of earthquakes should be extended in time, so that the seismic thermal anomalies should appear for at least two consecutive days. In that way, we extracted three thermal anomalies, which are Alarm I, II, and III, as shown in Figure 8. We then searched for the earthquakes that are potentially related to those thermal anomalies. According to the former cases [26,[34][35][36][37], we found that there is a strong spatial-temporal correlation between earthquakes with magnitude ≥5.5 and thermal anomalies extracted by the TFFA method. The latest thermal anomalies extracted by TFFA always occur within two weeks-i.e., a cycle of tidal force potential-before earthquakes, and the spatial distances between the thermal anomalies and earthquakes are smaller than or equal to 2 • . With these empirical conclusions, we selected earthquakes that occurred within 25 • -30 • N, 101 • -106 • E with magnitude ≥5.5 from December 2007 to December 2014 for this test. There are five such earthquakes, shown in Table 2 and as EQ1-EQ5 in the Figure 8.  In this study, the performance of thermal anomalies used as precursors is compared with random guessing. Random guessing means that there is no relation between earthquakes and thermal anomalies, and that they will occur at any time with equal probability. First, we get the prior probability that thermal anomalies at any time will correspond to earthquakes, assuming that the thermal anomalies are independent of earthquakes. For determining the correspondence between earthquakes and thermal anomalies, we set a temporal window of k days. For example, if the alarm A appears within k days before EQ1, the alarm A is judged to be related to EQ1. In return, using this temporal window and the earthquake catalog, we can determine a group of days such that, if these thermal anomalies occur (even randomly) on these days, they will correspond to earthquakes. For earthquake i occurring at t i , if alarms occur at times {t i − k, t i − k + 1, . . . , t i − 1, t i }, they will correspond to earthquake i and these alarms are regarded as true positives. Similarly, we can obtain N groups of days for all the N earthquakes, and Φ is the union of these N groups of days. We can then get the prior probability P for the alarms as follows: The number of true positives for alarms follows a binomial distribution: where h is the number of hit events or the number of true positives, and N A is the total number of alarms. Finally, we can get the P(h ≤ H ≤ N A ) as follows: Our null hypothesis is that thermal anomalies are not related to earthquakes. If P(h ≤ H ≤ N A ) < α, where α is the significance level, and if α << 1, we reject the null hypothesis, and thermal anomalies are related to earthquakes. In this study, we set α = 0.05. For the temporal window k, as is shown in Figure 4a, the earliest thermal anomalies regarded to be related to the Ludian earthquake appeared on 24 June, i.e., 40 days before the event. Therefore, we set k = 40 and P-value = 0.0112 < 0.05. We then reject the null hypothesis, and the long-term thermal anomalies around the epicenter of Ludian earthquake extracted by TFFA are related to earthquakes.

Discussion
Unlike the thermal infrared remote sensing data, the air temperature data provided by NCEP are not affected by clouds. When the area is obscured by clouds, it is impossible to extract thermal infrared anomalies caused by earthquakes, because thermal infrared data record cloud top radiation rather than land surface radiation. Moreover, the statistical methods for extracting thermal anomalies requires a lot of data processing, and the lack of good data will lead to bad and unreliable results. Figure 9 shows the cloud map provided by FY-2. The results show that the epicenter and even the whole study area were covered by clouds when there was a thermal anomaly under the clouds. In other words, it would have been impossible to get reliable, persistent, stable, and obvious thermal anomalies in the case of the Ludian earthquake using thermal infrared remote sensing data. It is concluded that, compared to thermal infrared data, air temperature data provided by NCEP have a huge advantage for monitoring the pre-earthquake temperature variation in cloudy areas. Strong cooling weather such as heavy rainfall and blizzard [41], however, can also hide the thermal anomalies, and how to extract the thermal anomalies before earthquakes under such weather needs further study. There is another huge limitation for the TFFA, because the key point of this method is the variation of tidal force potential at the hypocenter, but the information about the hypocenter is known only after the earthquake occurs. Therefore, this method is more suitable for studying the relationship between thermal anomalies and earthquakes, but cannot be used for earthquake prediction.

Conclusions
During cycles A, B, and C, all of the thermal anomalies are located on or near epicenters, and they generate from the bottom, spread upward, and disappear in the top layers. The vertical characteristic of those thermal anomalies is consistent with the results of the Mw 5.7 Jiujiang earthquake, and coincide with the principle of atmospheric thermal dynamic diffusion along the vertical: air heated by land, uplifted by heat flux, cooled and dissipated in the sky [26]. From period A to C, with the approaching earthquake, the thermal anomalies within each cycle become significantly more intense; moreover, the variation of the thermal anomaly within each cycle coincides with that of the tidal force potential. The thermal anomalies begin to rise at the quickly declining phase, reach a peak at troughs, and attenuate at the recovery phase. Based on the high spatial-temporal correlation between thermal anomalies and the Ludian earthquake, high consistency between spatial-temporal variation of thermal anomalies and tidal force potential, and the strong evidence provided by the significance test, it can be concluded that the thermal anomalies in periods A, B, and C are caused by the Ludian earthquake.
As shown in Figure 2, the variation of the tidal force potential is quasiperiodic. However, no earthquakes occurred during the similar phases at other times. Based on the opinion that the variance of air temperature can reflect the change of state of seismotectonic stress [26], we can propose that the relationship between the state of seismotectonic stress of cycles A, B, and C is as follows: C > B > A. In spite of the fact that there is energy accumulation caused by the tectonic displacement during periods A and B, it is released rapidly as the tidal force potential recovers. In period C, there are two rounds of rising air temperature. That means that the state of seismotectonic stress in period C increases until the critical breaking point of rocks, and then the earthquake occurs. The increasing seismotectonic stress during each cycle and the fact that no earthquake occurs during the similar phase of tidal force potential indicates that the state of seismotectonic stress is the decisive factor for triggering earthquakes, not the tidal force potential [26]. The Ludian earthquake is a good example to show the relationship between an earthquake, the tidal force potential, thermal anomalies (air temperature), and the state of seismotectonic stress; the tidal force potential can change the state of seismotectonic stress, which can change the local air temperature. However, earthquakes occur only when the seismotectonic stress reaches the critical breaking point of rocks. However, under the same tidal force potential, it is unknown what caused the different stress states in different periods, and what made the energy accumulation in period C much stronger than periods A and B. Exploring the potential factor is of great significance to the study of the mechanism of the Ludian earthquake and earthquake prediction.