Geothermal Anomalies and Coupling with the Ionosphere before the 2020 Jiashi M_s6.4 Earthquake

Abstract

Rock temperature reflects the adjustment in crustal stress, and the fluctuation of ionospheric electron concentration is closely related to short-term disturbances of the stress field. Their coupling may reveal short-term effects before strong earthquakes. This study explores the rock temperature changes and mechanical-electrical coupling in the lithosphere–ionosphere before the Jiashi M_s6.4 earthquake on 19 January 2020. The observed data were detrended by general polynomial piecewise fitting; three observation points within 150 km of the epicenter were found to show significant temperature fluctuations in the 15 days before the earthquake. The peak occurred synchronously five days before the earthquake, and the variation range was approximately 10⁻³ orders of magnitude. Five days before the earthquake, the electromagnetic satellite Zhangheng-1 synchronously observed an anomalous electron concentration in orbit near the epicenter, with a maximum value of 2.01 × 10¹⁰ m⁻³. The loading/unloading response ratio (LURR) was calculated using small earthquakes within 100 km of the epicenter; it showed that the large changes in rock temperature and the ionosphere occurred at high LURR, indicating high-stress accumulation in the region. Various anomalies appeared simultaneously and may indicate fault rupture, which may be caused by an acoustic-gravity wave, indicating a synchronous coupling between the lithosphere atmosphere and the ionosphere.

1. Introduction

Earthquake prediction is an important topic in earth science [1] and requires a deep understanding of the development of earthquake nucleation. The earthquake nucleation process has been extensively studied using various methods [2,3,4,5,6,7,8] to achieve breakthroughs in earthquake prediction. Many studies have reported that the tectonic stress field and its dynamic changes are key to the nucleation process of earthquakes [9,10,11]. The deformation and rupture of the lithospheric medium under load may change the underground temperature of the regional bedrock [12]. The physical theories of extrusion heating and tension cooling during stress loading are verified by rock mechanics experiments; they demonstrated that the tectonic stress variation could be identified from changes in bedrock temperature [13,14,15]. In practical earthquake research, Zoran [16] found that surface thermal anomalies developed before the 2011 M_s9.0 giant earthquake in Japan, and Qiang et al. [17] also detected significant thermal infrared temperature anomalies before the 1989 M_s6.1 earthquake in Datong, Shanxi, China. Green [18] and Yasuyu et al. [19] reported that sliding friction caused by the rapid dislocation of the fault plane during the earthquake was the main reason for the sharp increase in temperature near the fault plane. Chen et al. [20] studied the temporal and spatial evolution of the latent heat flux on the Earth’s surface before the Chile M_s8.8 earthquake and also discussed the theoretical relationship between latent heat flux and surface temperature. The research shows that the strain, friction, and fracture of rock, which are caused by the change of stress, may change the temperature. The observed temperature change may help us understand the stress adjustment before an earthquake.

Ground temperature changes have typically been observed by satellite thermal infrared remote sensing and other imaging technologies. Ma et al. [21] studied the possibility of using thermal infrared anomalies to detect incoming earthquakes on the basis of the relationship between infrared brightness temperature and terrestrial heat flow. Chen et al. [22] used surface temperature data from the Terra and Aqua satellites to study the co-seismic thermal response of the Wenchuan M_s8.0 earthquake. Choudhury et al. [23] and Guo [24] studied pre-earthquake thermal anomalies using remote sensing technology. These studies showed that thermal infrared data on geothermal changes might contain information on the occurrence of earthquakes. However, observation accuracy and excessively uncertain space range are key scientific issues that urgently require a solution. The altitude of satellite thermal infrared observation is generally more than 30,000 km, with the anomaly covering a range of hundreds to thousands of kilometers [23] and the observation accuracy of the thermal imager reaching 0.01 degrees. However, the altitude of the ionospheric satellite observation is just 500 km, with an anomaly coverage area of several hundred kilometers. Thus, the ionospheric anomaly is selected to couple with the rock temperature observation to explore earthquake preparation [25,26,27].

Recent studies have suggested that earthquakes are closely related to bedrock fractures [28], and arranging a certain number of temperature sensors in key areas may be an effective way to detect temperature changes in bedrock before earthquakes. In fact, Chen et al. [29,30,31] demonstrated the feasibility of bedrock temperature observation of earthquakes and analyzed changes in rock temperature observed before the Lushan M_s7.0 and Kangding M_s6.3 earthquakes. The pre-earthquake variation in rock temperature is thought to be related to variation in the tectonic stress field, which reflects regional stress enhancement. Using wavelet transform analysis, Zeng et al. [32] identified a rapid fluctuation change in ground temperature of 0.001–0.002 °C approximately 80 days before the Lushan M_s7.0 earthquake. Changes in bedrock temperature indicate the state of the rock during stress adjustment, which is part of earthquake formation. These precursory anomalies in ground temperature provide a reference for understanding the development of strong earthquakes.

Earth is a dynamic system, and changes in the lithosphere will probably cause atmospheric resonance. The coupling of multiple layers before strong earthquakes has attracted attention. Chen et al. [33] integrated observational methods for atmospheric pressure, vertical electric field, atmospheric radon concentration, groundwater level, and precipitation to construct the MVP-LAI (monitoring vibrations and perturbations in the lithosphere, atmosphere, and ionosphere) system, which they used to observe the surface vibration, atmospheric pressure, and electron concentration before the 2021 Luxian M_s6.0 earthquake synchronous coupling effect. Using ground electromagnetic data, Li et al. [34] analyzed the electromagnetic coupling effect three days before the 2008 Wenchuan M_s8.0 earthquake. Carbone et al. [35] constructed a coupling model for the period before strong earthquakes using ground motion and ionospheric data. Notable changes in ionospheric data before strong earthquakes have been extensively studied. They were likely caused by the acoustic-gravity wave or electric field generated in the lithosphere before the earthquake [33,36], but further research may be needed [37]. The communication of energy and matter exists between the lithosphere, atmosphere, and ionosphere. If the coupling relationship between multiple layers can be further understood, it will enrich our study of earthquakes.

Before the M_s6.4 Jiashi earthquake in Xinjiang in 2020, there was an anomaly in bedrock temperature [38,39], caused by changes in crustal stress and secondary fluid effects. In this work, using observational data on rock temperature and the ionosphere, the relationship between the various characteristics of rock temperature before earthquakes and the sub-instability process are studied using cross-layer coupling. It is found that the stress release of the fault area accelerates and changes from quasi-static to dynamic, which is similar to the variation of faults before instability studied by Ma et al. [40].

2. Regional Tectonic Background

On 19 January 2020, the Jiashi M_s6.4 earthquake occurred in the Kepingtag thrust fault. The intense fault activity and the frequency of large earthquakes cause the fracture of the alluvial fan, forming a series of fault scarps at the foot of the mountain. Kepingtag tectonic thrust belt is in the western segment of the southern Tianshan Mountains. The geothermal resources in the South Tianshan Mountains are abundant, and the geothermal distribution is closely related to the tectonic system. The geothermal gradient in the center of the anomaly is 2–3 times higher than the normal gradient. The terrestrial heat flux in the Southern Tianshan Mountains is about 70–75 mW/m², and the adjacent Pamir region is about 90–150 mW/m² [41]. The belt is located in the middle section of the Indian to Eurasian plates and is among the youngest active intracontinental orogenic belts [42]. The Global Navigation Satellite System (GNSS) velocity field shows that the movement rate of each station maintained the same direction (Figure 1), and the average near-north–south velocity was approximately 17.6 mm/a. The velocity differs significantly between the north and south sides of the fault [43], indicating relatively high-stress accumulation in this area. The compressive stress generated by the collision of the Indian and Eurasian plates is transmitted over a long distance, and the lithosphere of the Tarim Basin is subducted below the Tianshan Mountains, resulting in the uplift of the mountains. The large crustal stress accumulation caused by subduction leads to frequent strong earthquakes in this area [44,45,46,47,48]. Since 1990, 15 earthquakes with M_s ≥ 6.0 (from China Earthquake Networks Center) have occurred within 200 km of the epicenter of the Jiashi earthquake.

3. Observational Data Sources and Processing

3.1. Source and Processing Method of Rock Geothermal Observation Data

From 2016 to 2019, the Institute of Geology of the China Earthquake Administration and the Earthquake Administration of the Xinjiang Uygur Autonomous Region jointly established a number of bedrock temperature monitoring stations along the Keping fault (F1 in Figure 1) to detect regional fault stress changes and earthquake development processes. The bedrock temperature is monitored using wireless geothermal telemetry equipment developed by the Institute of Geology, China Earthquake Administration, with a sampling period of 10 min and an accuracy of 0.00003 °C. Holes are drilled in the bedrock, and a temperature sensor is installed in each hole using high-grade cement, which is closely fitted to the bedrock (Figure 2). The sensors are vertically distributed at multiple levels between the surface and 50 m underground. The observation stations are selected along the fault, with relatively complete bedrock and far away from water sources. The topographic environment is similar, without special terrain. The observed bedrock temperature data are stable and smooth, and the layer below 10m is not directly affected by the air temperature. The borehole is sealed by a cement tank, which is hard to be affected by external interference factors. The observed data reflect the seasonal variation of shallow ground temperature, but the deep ground temperature is constant.

In this study, the bedrock temperature observation data were processed by general polynomial piecewise fitting. To ensure the continuous smoothness of the fitting curves of adjacent segments, the following constraints were introduced: the r-order derivatives of the fitting curves in adjacent data segments at the common points are equal, r = 0, 1, 2..., and r is less than the order of the fitting polynomial. The data are divided into p data segments, and each segment is fitted by m polynomials. The total number of parameters to be calculated is thus p × (m + 1); the constraint condition is (p − 1) × m, and the independent parameters are

$$r=p\cdot (m+1)-(p-1)\cdot m=p+m$$

(1)

The required solution is obtained by polynomial fitting of the r parameters. Because the continuous smoothness of the curve at the point should be considered, the required curve is written as

$$y(x)={\displaystyle \sum _{i=0}^m{a}_i{x}^i}+{\displaystyle \sum _{j=1}^{p-1}{a}_{m+j}{(x-{\xi }_j)}_{(+)}^m}$$

(2)

$${(x-{\xi }_j)}_{(+)}^m=\left\{\begin{array}{cc}{(x-{\xi }_j)}^m& \begin{array}{cc}& (x\ge {\xi }_j)\end{array}\\ 0& \begin{array}{cc}& (x<{\xi }_j)\end{array}\end{array}\right.$$

(3)

where ${\xi }_j$ is the x value at the j-th knot point. From formula (2), we obtain

$$\left.\begin{array}{l}{f}_{j+1}=f_j(x)+a_{m+j}{(x-{\xi }_j)}^m\hfill \\ {f^{\prime }}_{j+1}(x)={f^{\prime }}_j(x)+m{a}_{m+j}{(x-{\xi }_j)}^{m-1}\hfill \\ {f^{″}}_{j+1}(x)=f^{″}(x)+m(m-1)\cdot a_{m+j}{(x-{\xi }_j)}^{m-2}\hfill \end{array}\right\}$$

(4)

When $x={\xi }_j$, the second term of the above equations is 0, and the function value of the fitting curve at the point is equal to the first and second derivatives. That is, the fitting curve is continuous and smooth at each point $(x={\xi }_j)$. We treat $x-{\xi }_1,x-{\xi }_2,\cdots ,x-{\xi }_{p-1}$ as a new observed quantity and take $a_{m+1},a_{m+2},\cdots ,a_{m+p-1}$ as the new undetermined coefficient. Least-squares processing of $p+m$ undetermined coefficients is performed to solve for all fitted observed data samples n. The residual error equation is

$$v_k=y_k-{\displaystyle \sum _{i=0}^m{a}_i{x}_k^i-{\displaystyle \sum _{i=1}^{p-1}{a}_{m+j}{(x_k-{\xi }_j)}^m}}$$

(5)

The least-squares method is used to obtain all the fitting coefficients according to the least-square sum of errors. The regular equation can be solved by the Gaussian reduction method. Finally, all fitting coefficients, fitting values of the curve, and fitting residuals and their standard deviations are the output [49]. This method eliminates the effect of erroneous data and yields information on changes that occur before earthquakes.

3.2. Ionospheric Data Sources and Processing Methods

Electron density data from the Langmuir Probe (LAP) of the electromagnetic satellite Zhangheng-1(known as CSES) are used as ionospheric data. The satellite was launched on 12 February 2018. It is a sun-synchronous satellite with an orbital inclination of 97.4°, a flight altitude of 507 km, and a revisit period of five days. The indication of local time is about 2 AM for Zhangheng-1. The LAP makes measurements every 3 s in survey mode. When passing over the seismic zone and Chinese region, it switches to detailed inspection mode and makes measurements every 1.5 s. Because the ionospheric changes at night are less affected by interference, this study analyzed only the ascending orbit data. These data correspond to local night, and data with Kp ≥ 3 (Kp is an index that describes global geomagnetic activity) are excluded to eliminate the effects of magnetic disturbances. According to previous studies, most ionospheric disturbance anomalies occurred within 15 days before an earthquake [50,51]. The ascending orbit data for 30 days before the earthquake were selected for analysis, and the Kriging method (a method for unbiased optimal estimation of variables in a finite region based on variogram theory and structural analysis) was used to interpolate between adjacent orbits.

Due to the ionospheric electron density exhibiting relatively significant annual and seasonal changes, the ascending orbit data in the quarter before the earthquake were used as background data, where data with Kp ≥ 3 were removed. A 5° × 5° grid was used, and the average value of each grid in the quarter was calculated as the background value to analyze changes during the study period in this area. According to the orbital position of the electromagnetic satellite, the orbit closest to the earthquake was found, and the part that was closer to the epicenter and had higher abnormal values was selected for the study.

4. Changes in Rock Temperature before the Earthquake

4.1. Fluctuations Immediately Preceding the Earthquake

The raw data on bedrock temperature is relatively stable. After the data were detrended using MAPSIS software, the three observation points at Gedaliang, Siker, and Xiaohaizi began to show clear synchronous fluctuations with the same shape 15 days before the earthquake. Note that the sensor at the Machang observatory is different from others, with an accuracy of 0.0005 °C, which is lower than that of the Gedaling, Siker, Xiaohaizim, and Wupaer stations (the accuracy is 0.00003 °C). Thus, the quality of the observed curve at Machang is lower than the rest of the stations (Figure 3). Due to the large distance and different tectonics, no similar pre-seismic changes were observed in the Machang/Wupaer. The maximum fluctuation occurred on 14 January 2020, with an amplitude of approximately 0.005 °C (the green shadow zones of Figure 3). According to the sub-instability theory [40,52,53,54], this regional bedrock temperature change may indicate that the fault was entering sub-instability. The change of five days before the earthquake may indicate information about the earthquake. This time is almost the same as the perturbation in water flux recorded by De Luca et al. before the 2016 Amatrice (Italy) M_w6.0 earthquake [55]. The temperature sensor is sealed in the underground rock and is hard to be affected by external factors. Moreover, each sensor is independent and cannot affect the others. Other influences, such as voltage fluctuation, instrument failure, and subsurface fluid migration, can also be excluded because the voltage and instrument are normal during the temperature change, and the observed data are stable and smooth.

4.2. Short-Term Changes

To explore the correlation between bedrock temperature changes and earthquakes, we used the raw data dynamic tracking method to analyze the observed bedrock temperature data within 200 km of the epicenter. Relevant temperature anomaly effects occur within 200 km around the epicenter before M_s≧ 6 earthquake [56]. The observation results are shown in the green shadow zones of Figure 4. These data changes relative to the daily data are more prominent. The changes essentially reflect the bedrock temperature change. These changes began 120 days before the earthquake.

The Wupaer and Aqia observatories measured rapid increases in September 2019. The Machang and Siker observatories showed rapid increases in October 2019. The Gedaliang observatory experienced a rapid increase in November 2019. These changes only occur in some layers of the observation station, which may be related to the lithology and core integrity of each layer. No similar changes were found in the Xiaohaizi observation station. The time, amplitude, and morphology of the short-term temperature change look different from that of the pre-seismic fluctuations, which may be caused by stress change, underground fluid, and other factors. The pre-seismic fluctuations occurred almost simultaneously, while the time of the short-term changes was quite different. The range of the pre-seismic fluctuations is about 10⁻³, while the range of the short-term variation is 10⁻¹. In addition, the trend turnings, significant increases or decreases in the short-term changes are different from the pre-seismic changes. Figure 5 shows the evolution of the epicentral distance of the temperature increases observation stations with time, which clearly shows a gradual shift toward the epicenter. This decrease in the epicentral distance further suggests the correlation between the ground temperature change and the seismogenic process.

5. Ionospheric Observation Anomalies

In agreement with the timing of the temperature fluctuations before the earthquake (Figure 3), the electron density measured by the electromagnetic satellite Zhangheng-1 also showed an abnormally high value on 14 January 2020, with the anomaly region indicated by Letter “A” (Figure 6b). In the ascending orbit, the satellite observed the region closer to the epicenter and measured abnormally high values in the study area, which has latitude and longitude ranges of 79.0–81.5° E and 38.0–46.0° N (the thicker part of the orbit in Figure 6b). The electron density time series in Figure 6a indicates relatively stable electron density on this orbit from 5 December 2019 to 9 January 2020, with fluctuations of approximately 1.2 × 10¹⁰ m⁻³. However, on January 14, the maximum average value was 2.01 × 10¹⁰ m⁻³, which exceeded the background value by 67%; this observation was synchronous with the rock temperature anomaly. The variation of spatial electron concentration is generally related to solar activities and near-surface disturbances, such as cyclones and volcanic eruptions. Earthquakes can also affect it through atmospheric fluctuations, electric field disturbances, and magnetic field disturbances. Before the earthquake, due to stress accumulation, rock fracture, fissure increase, gas escape, greenhouse gas release, heat source, and others, the near-surface atmosphere, electric field, and magnetic field should be changed. The electron concentration oscillates when it reaches the ionosphere, which may be caused by acoustic gravity waves [35,57].

6. Loading/Unloading Response Ratio Analysis

He et al. [58] defined the stages of earthquake nucleation as stress accumulation, initial weakening, and overall weakening, where core weakening appears in the overall weakening stage. The boundary of the core weakening zone initially expands with time and then contracts in the later period. Finally, it points to the fault block, which is consistent with the changes in rock temperature before the earthquake. It also confirms the transformation of the sub-instability stage from the independent activity of each part of the fault to the overall activity, which was defined by Ma et al. [52] on the basis of rock experiments. The overall activity of the fault results from high stress. To verify the high stress at this time, we used the loading/unloading response ratio (LURR) method. The LURR quantitatively reflects the approaching instability of a nonlinear system [59]. When the regional stress is in the normal state, the earthquakes in the loading and unloading stages should be about the same. When the region is in a high-stress state, the stress disturbance will cause an earthquake more easily. It is calculated as follows:

$$Y=\frac{X_{+}}{X_{-}}$$

(6)

where the plus sign (+) and minus sign (−) indicate loading and unloading, respectively. The response rate X is expressed as

$$X=\underset{\Delta P\to 0}{\lim }\left(\frac{\Delta R}{\Delta P}\right)$$

(7)

where Δ R and Δ P represent the changes in the load P and response R under a load P. The formula for using seismic energy E as a response is as follows:

$$Y=\frac{{\left({{\displaystyle \sum }}_{i=1}^{N_{+}}{E}_i^m\right)}_{+}}{{\left({{\displaystyle \sum }}_{i=1}^{N_{-}}{E}_i^m\right)}_{-}}$$

(8)

When m = 1, E^m represents the seismic energy, which can be calculated according to the Gutenberg formula; for m = 1/2, E^m represents the Benioff strain. In this study, we use the Benioff strain to calculate LURR. N⁺ and N⁻ represent the number of earthquakes during loading and unloading, respectively. We classify the earthquakes within 100 km of the epicenter by loading and unloading (Figure 7a) for 200 days before the earthquake. Then we selected the 0 < M_s < 4 earthquake catalog to calculate LURR. The focal mechanism parameters were selected from the moment tensor results calculated using the Global CMT page: strike: 80°, dip: 71°, slip angle: 124°, focal depth: 12 km. The friction coefficient of the fault was set to 0.4. A window length of 30 days and a step size of 10 days were used to calculate the LURR. The LURR was found to be abnormally high on 31 December 2019 (Figure 7b), and the Jiashi M_s6.4 earthquake occurred on 19 January 2020. A high-value anomaly also occurred on 8 October 2019; this anomaly preceded the Wushi M_s5.0 earthquake on 27 October 2019 by 19 days. The magnitude increased with increasing LURR. Figure 7b shows that the ratio is typically approximately one and reaches six before the Jiashi earthquake, indicating a state of high stress before the earthquake. Which indicated that the regional stress level was higher and the load triggered more earthquakes. The green area in Figure 7a shows that the load earthquakes increased.

7. Discussion

From 120 to 60 days before the earthquake, the bedrock showed different degrees of abnormally high temperature, which shifted dynamically, appearing closer to the epicenter over time, indicating that the ground temperature changes reflected the development of the seismogenic region. Zhang et al. [60] found a similar anomaly in geophysical observation data before the Jiuzhaigou M_s7.0 earthquake in China in 2017. Shimojo et al. [61] also found that mild earthquake activity before the earthquake accelerated toward the epicenter in a study of the nucleation process of the Changye M_s6.7 earthquake in Japan in 2011. The regional anomaly beginning far from the epicenter and approaching it may be generated by crustal stress, reflecting the action of the regional stress field, which is also indicated by the 3.0 ≦ M_s≦ 3.9 apparent stress in the Keping block in the same period [62]. This indicates that the crustal stress adjustment information can be reflected through the change in ground temperature, and the crustal stress adjustment reflects the spatial-temporal process of earthquake preparation, which can be applied to detect the occurrence of earthquakes.

The changes in bedrock temperature appeared at three observation points (Gedaliang, Siker, and Xiaohaizi) five days before the earthquake, and the peak fluctuation appeared synchronously. The fluctuations began approximately 15 days before the earthquake, and the amplitude was accelerated. The rock temperature both increased and fluctuated before the earthquake. This change in temperature may indicate that the rock mass is also in a state of acceleration and vibration. According to the studies presented by Li et al. (2020), there were also marked anomalies in the gravity high-gradient zone before the earthquake [63]. It is speculated that the fault had entered a stage of acceleration and vibration at this time, which is the key to sub-instability [40]. Sub-instability indicates increased stress. The LURR reflects dynamic changes in the meso-constitutive relation of the rock in the source area. When the rock system is in the elastic stage, the LURR is low (approximately 1); when the system is in the expansion stage, the LURR is greater than 1. When the peak loading/unloading response is reached, the rock system may have entered a sub-instability stage. At this time, the LURR decreases. As shown in Figure 7b, the peak LURR appeared 19 days before the earthquake, and the synchronous coupling between rock temperature and the ionosphere five days before the earthquake may indicate that the rock shifted from the steady state to the sub-instability stage (BC segment in Figure 8). The rock change rate is accelerated by the slope of the straight line L2 in Figure 8. The constitutive curve of rock shows the change of rock during stress strengthening. When the stress is low, the system is in the elastic stage, and the changes are recoverable. When the stress is high, the rock is in the dilatant process. With the development of cracks, the deformation is irrecoverable. When the stress reaches point B, it exceeds the limit of the rock, and the system enters the stage of sub-instability. In this stage, accelerated generation of cracks can be observed. In the sub-instability stage, the micro-vibration in the lithospheric medium causes the rock temperature to shift, and the change in the surface heat source and escaping gas released form acoustic-gravity waves in the atmosphere; these waves enhance gravity waves, which propagate to the ionosphere and cause synchronous oscillatory changes in electron concentration [57]. The change of 5 days before the earthquake is consistent with the abnormal change time of Yan Rui et al. based on ionosphere observation [64]. According to this theory, the change in rock temperature in the area before the Jiashi M_s6.4 earthquake and the accelerated vibration of the fault may have caused the change in the thermal field of the surface atmosphere. The raised temperature in rock accelerates the release of underground chemical gases, including greenhouse gases. Second, heating the near-surface atmosphere causes atmospheric fluctuation. These may cause a local electric field change or the formation of an acoustic gravity wave, resulting in a change in ionospheric electron concentration. The propagation of acoustic gravity waves is affected by Coriolis force, wind speed, nonlinear propagation, and other factors. The ionospheric electron concentration anomaly distribution deviates from the epicenter area. However, the results presented in this paper show that the ionospheric electron concentration anomaly and the changes in the acceleration of rock vibration before fracture are both associated with the preparation for large earthquakes.

8. Conclusions

This study showed that information on seismic nucleation and rupture could be captured by observing changes in rock temperature and their coupling with the ionosphere. Observations of the Jiashi M_s6.4 earthquake showed that the bedrock temperature is sensitive to changes in stress loading and exhibited space-time consistency. Mechanically, the seismogenic characteristics are consistent with the transition from the elastic stage to the plastic yield stage. The spatial-temporal evolution of the medium-term changes of rock temperature can be applied to determine the location of future large earthquakes, while the short-term synchronous changes may reveal the short impending earthquake potential in the region. The obvious acceleration and fluctuation of ground temperature at multiple observation stations before the earthquake demonstrated that the fault was entering the sub-steady state. The time series of the LURR confirms the relationship between ground temperature variation and fault sub-instability. This accelerated fluctuation should be generated in the context of strong stress. In addition, the vibration of the fault and the change in the temperature field may form an acoustic gravity wave, causing the ionosphere to oscillate synchronously. The synchronous changes in rock temperature and the ionosphere may reveal fault instability. These findings are useful for further exploration of the coupling mechanism of the geophysical field before strong earthquakes and the development mechanism of earthquake nucleation and rupture. However, the coupling mechanism of rock temperature and the ionosphere requires further verification by theoretical models.