Multi-Parametric Climatological Analysis Reveals the Involvement of Fluids in the Preparation Phase of the 2008 Ms 8.0 Wenchuan and 2013 Ms 7.0 Lushan Earthquakes

Abstract

A multi-parametric approach was applied to climatological data before the Ms 8.0 2008 Wenchuan and Ms 7.0 2013 Lushan earthquakes (EQs) in order to detect anomalous changes associated to the preparing phase of those large seismic events. A climatological analysis for seismic Precursor Identification (CAPRI) algorithm was used for the detection of anomalies in the time series of four parameters (aerosol optical depth, AOD; skin temperature, SKT; surface latent heat flux, SLHF and total column water vapour, TCWV). Our results show a chain of processes occurred within two months before the EQs: AOD anomalous response is the earliest, followed by SKT, TCWV and SLHF in the EQs. A close spatial relation between the seismogenic Longmenshan fault (LMSF) zone and the extent of the detected anomalies indicates that some changes occurred within the faults before the EQs. The similarity of time sequence of the anomalies between the four parameters may be related to the same process: we interpret the observed anomalies as the consequence of the upraising of gases from a fluid-rich middle/upper crust along pre-existing seismogenic faults, and of their release into the atmosphere. Our multi-parametric analytical approach is able to capture phenomena related to the preparation phase of strong EQs.

1. Introduction

Earthquake preparation is a long-term lithospheric process and it is driven by tectonic- and/or fluid-induced stress changes [1,2,3,4]. The earthquake preparation processes are, however, poorly constrained and far from being fully understood [5,6,7]. To date, all precursor information has been observed on the surface, near surface, atmosphere, and ionosphere away from the hypocenter of earthquakes (EQs). How can this information be transferred from the depth of focus to the location of the observing instrument? What kind of changes may occur in the transfer process? What interactions exist between different physical and chemical information? It is fundamental to deeply understand the process of EQ occurrence and explore new scientific ideas of EQ prediction. Based on actual observation data combined with certain scientific assumptions, it is the necessary measure to advance the scientific development of EQ prediction, by carrying out multidisciplinary comprehensive exploration, identifying multiple physical and chemical phenomena, and their transfer and diffusion processes.

A lot of effort has been put into analyzing a wide variety of phenomena preceding seismic events with the aim of finding some recurrent and recognizable patterns: among others, groundwater level changes, gas and infrared (IR) electromagnetic emissions, atmospheric optical depth (AOD) changes, thermal anomaly, surface deformations, electric and magnetic fields [5,8,9]. As an example, thermal anomalies on seismogenic faults are sometimes observed starting from several months to a few days before strong earthquakes (EQs). The changes observed during the preparation phase of large seismic events are generally explained by the lithosphere–atmosphere–ionosphere coupling (LAIC) model [7], which suggests that a diversity of phenomena might occur before EQs due to either an exchange of energy or particles between the lithosphere and the atmosphere reaching, in some cases, the ionosphere. As an example, heat fluid could be split into water and gas, at depths of a few kilometres. Air temperature increases caused by the latent heat release as a result of water vapor condensation on ions, these latter produced by air ionization due to radon. The rapid lithospheric degassing would trigger several atmospheric processes near Earth’s surface, namely: air molecule ionization and plasma chemical reactions, which dominate the formation of complex molecular ions and of large ion clusters up to aerosol size. Traditionally, meteorology considers aerosols as the main centers of nucleation and water condensation [7]. Therefore, the heat, water vapor and other gas could rise to the surface of the Earth and then the lower atmosphere through the geological structures (as faults, cracks, fractures etc.), thus changing other parameters, such as AOD. Therefore, considering geochemical–thermal anomaly interaction of the LAIC model as a theory, the objective of this study was to observe the abnormal response of the multi-parameter data from the lithosphere and atmosphere in the seismogenic process and to further explain the EQ mechanism along the Longmenshan fault (LMSF) zone. In this study, we selected four climatological parameters: skin temperature (SKT), total column water vapor (TCWV), surface latent heat flux (SLHF) and aerosol optical depth (AOD) to represent the geochemical–thermal process of the lithosphere and atmosphere prior to EQs.

Regarding the anomaly phenomena, thermal infrared radiation (TIR) anomalies were observed by satellites before EQs [10]. The increase of skin temperature before the EQ was also regarded as a possible anomaly due to tectonic activities [11]. Other thermal parameters, such as the outgoing long-wave radiation (OLR) and surface latent heat flux (SLHF), have shown anomalous changes [12,13]. The latter (SLHF) can be explained as the gas emission increases the heat flux convection and may lead to the increase of SLHF. Ganguly (2016) reported an increment higher by 40% and 6% of AOD and columnar ozone, respectively, before Nepal EQ on 25 April 2015 [14]. In association with the strong (Mw7.8) EQ that occurred in Ecuador on 16 April 2016, variations of the AOD obtained from MODIS (MODerate-resolution Imaging Spectroradiometer) onboard the AQUA satellite, were observed to peak 10 days before the event [15]. AOD abnormal phenomenon in Wenchuan EQ was observed by Qin et al. [16]. Twelve days before the mainshock, a strip-shaped high-value AOD anomaly was also found at the north Iraq-Iran border [17]. Therefore, the analysis of climatic parameters could, in principle, give us information on the preparatory phase of large seismic events and on the EQ triggering factors.

Here, we present an analysis of climatological parameters in a long time period preceding the Ms 8.0, May 21 2008 Wenchuan and Ms 7.0, 20 April 2013 Lushan EQs. The two EQs were the most significant events in the region in the past 12 years and have great regional representativeness. These events occurred in the eastern margin of the Tibetan Plateau, in the LMSF zone, at 20 and 18 km depth, respectively (Figure 1). The 2008 Wenchuan EQ produced a 240-km surface rupture along the Beichuan-Yingxiu Fault (BYF), a 70-km rupture along the range front Guanxian-Jiangyou Fault (GJF), and a 6-km rupture along the NW-trending arcuate Xiaoyudong Fault [18,19,20]. The distribution of Wenchuan EQ aftershocks after their relocation evidenced a total rupture length of about 330–350 km [21]. Many people were involved directly (more than 69,000 deaths, 374,176 injured) and indirectly (about 30 million people living in 16,000 square kilometers near the epicenter were affected).

Lushan EQ occurred in Lushan County, Sichuan Province, about 100 km away from Chengdu. Direct economic losses were estimated to be more than ten billion yuan (RMB). Its epicenter laid nearby the Shuangshi Dachuan fault (SDF), which is the southern segment of the LMSF zone. The focal depth of the mainshock was 13 km. The 2013 Lushan EQ occurred within a blind thrust fault located between the Shuangshi-Dachuan Fault and the Dayi Fault [22]. The relocations of the Lushan aftershock sequence indicated a rupture length of about 35 km. The earthquake killed nearly 200 people, caused more than 10,000 injuries, and damaged many buildings [23,24,25].

The methodology of this study consists of the integrated investigation of several atmospheric and environmental parameters in search of anomalous variations which have been found to be possibly related to approaching EQs [26,27]. We considered a fixed radius of 700 km around the epicenters and long-time series (about 30 years). Now, many studies on the thermal anomalies are based on the satellite data; however, the analysis is not an easy task because opaque clouds may mask surface temperature, gas, and aerosol optical depth data. Although there are a large number of ground station data and satellite monitoring data, they are not evenly distributed geographically. At the same time, some data have not observed all or some features in some areas, and even some data might be wrong. The methodology of “Reanalysis data” merges the observation data from various sources, different error information, and different spatial and temporal resolutions into the numerical dynamic model. According to strict mathematical theory during the reanalysing process, an optimal solution is found between the model solution and the actual observation to assure the accuracy of the data. Through the retrospective analysis of different climatological quantities, long-term continuous and spatially complete datasets with numerical assimilation models were produced. Compared with the complete satellite data, the “reanalysis data” have the great advantages of completeness, resolution and accuracy. Therefore, “reanalysis data” from European Center Medium Weather Forecast (ECMWF) of SKT, TCWV, and SLHF and the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2) of AOD were used in this study. Since it was reported that the occurrence of anomalies were noticeable 3 months before the Lushan EQ [28,29,30], for both Wenchuan and Lushan events, we analysed the time series of all parameters starting three months before the corresponding mainshock occurrence. By a multi-year meteorological data analysis, the greenhouse effect needs to be considered, which results in the temperature increase.

In the following sections, we describe the seismotectonic and geological setting, the used data, the analytical methods and the applied algorithms. In the last section, we discuss the results by integrating our results with the available geological, geophysical and geochemical information. We provide evidence that the upraising of crustal fluids played a significant role in the preparing phase of the 2008 Wenchuan and 2013 Lushan EQs and that the long-term climatological analysis of seismogenic zones may give information on the processes occurring in the lithosphere.

2. Tectonic Setting

Wenchuan and Lushan EQs occurred in the LMSF region, in the eastern margin of the Tibetan Plateau, which is uplifting and expanding towards East. The Wenchuan EQ was the first seismic event of such large magnitude recorded on a high-angle (60-80°) thrust fault within a continental interior [31]. The Lushan event occurred on a low angle reverse LMSF thrust between the Bayan Har Block and the South China Block [32]. Focal depth profiles of Lushan EQ show fault planes dip to northwest as a low-angle thrust [23]. NE–SW-striking fault-related folds and nappes associated to the thrusts include Precambrian basement rocks and accommodate most of the crustal shortening [24]. This regional scale geodynamics and the associated deformation make the Tibetan Plateau and its surrounding area one of the most active region of crustal deformation and seismicity in China [33]. The LMSF zone is a first-class boundary of active blocks and huge potential structure background for strong EQs and consists of different minor segments. This major fault system separates the Tibet chain from the Sichuan basin. The aftershocks of the Wenchuan and Lushan EQs propagate from 18-20 km to the surface along the Guanxian-Anxian Fault and Yingxiu-Beichuan Fault segments of the major LMSF zone [34]. A 30–40 km wide coseismic gap zone separates the the 2013 Lushan and 2008 Wenchuan EQs [35]. Tomographic images show that the Lushan and Wenchuan EQs occurred in two areas characterized by high velocity (Vp, Vs). These areas probably reflect the Precambrian meta-sediments and/or igneous rocks. The low velocity zone separating the hypocenters of the two EQs, which extends in depth from the middle-crust to the surface corresponds to a high conductivity zone revealed by magnetotelluric data [36]. Zhao et al. (2012) reported a low-resistivity zone in correspondence of the LMSF area at a depth of ~20 km beneath the eastern Tibetan Plateau that terminates ~25 km west of the responsible for the Wenchuan EQ [37]. According to Lei and Zhao (2009), the fluid-filled fractured rocks affecting this area and the LMSF zone may have influenced the generation of the large Wenchuan EQ [38]. Cui et al. (2017) reported that large amounts of crustal CH₄ and CO₂ were released into the atmosphere before, during and after the Lushan and Wenchuan EQs, along the fault zones [28]. These anomalies may be due to porosity changes in the fluid-filled fractures of the thrust before and after the seismic events.

In this study, the investigated areas of the Wenchuan and Lushan EQs fall within the range of six degrees from the epicenter.

3. Materials and Methods

3.1. Climatological Data

(1) ECMWF

The European Centre for Medium-Range Weather Forecasts (ECMWF) is a research institute which produces numerical weather predictions at global scale (source: http://www.ecmwf.int/en/about).

ECMWF uses its prediction models and data assimilation systems to "re-analyze" the observations produced by atmospheric models, creating global datasets that describe the recent behavior of the Geosphere (lithosphere, hydrosphere and the atmosphere). The aforementioned data were used to monitor for research, e.g., climate change and also for commercial use [39].

ERA-Interim is one of the projects in the global atmospheric reanalysis that starts from 1979 and is updated in real time. The data assimilation system that ERA-Interim produces uses a unique integrated model of ECMWF atmospheric forecasting for the whole time series (December 2006 version). Data assimilation is based on a 12-hour four-dimensional variational analysis, using an adaptive bias estimation for the satellite radiation data.

In this study, SKT, TCWV and SLHF ERA-Interim datasets were analyzed starting from 1980 to 2007 for Wenchuan EQ, or 2012 for Lushan EQ, and compared with the same period of 2008 or 2013, respectively, to observe the results.

As regards the year in question, it was decided to use the ECMWF analysis data, the operational archive, in order to avoid possible effects of excessive homogenization of the data caused by the ERA-Interim reanalysis procedure, which could have hidden the potential precursor anomalies with dataset biases.

Values at synoptical date, 00:00, 06:00, 12:00 and 18:00 UTC were collected. Since thermal anomalies related to normal weather variations are scarce at night-time, in the absence of solar radiation contribution, only local night-time data were used, i.e., at 18:00 UTC, corresponding to around 2 am LT. This permits to distinguish possible temperature anomaly variation induced by earthquake from the usual weather phenomena.

The spatial resolution used is 0.5° by 0.5° grid box, over a circular area of 700 km of radius.

(2) MERRA-2

“Modern-Era Retrospective analysis for Research and Applications, version 2” (MERRA-2) dataset was also used: MERRA-2 is a NASA atmospheric reanalysis for the modern satellite era (1980 – onward with a delay of one-two months), using the Goddard Earth Observing System Model, Version 5 Atmospheric Data Assimilation System (ADAS), version 5.12.4.

MERRA-2 is also the first satellite-era global reanalysis to assimilate space-based observations of aerosols, including assimilation of bias-corrected AOD from AVHRR and MODIS, MISR AOD over bright surfaces, and AERONET AOD, by using the Goddard Aerosol Assimilation System (GAAS) [40].

In this study, M2I3NXGAS, a MERRA-2 subset, was used. It is an instantaneous, single-Level (vertically integrated), chemical and aerosol model of atmosphere with 3-hourly values supplied in NetCDF format. Spatial resolution is 0.5° by 0.625° grid box, in latitude and longitude, respectively.

3.2. Analysis Method

Based on the reanalysis data used in this study, the time series of analysis data were up to 33 years. Taking the parameter of SKT as an example, the 33-year time series will inevitably have a rising trend brought by “global warming”. The CAPRI method is used mainly to remove a possible ‘‘global warming’’ effect, avoiding to classify as abnormal a more recent year just because of global warming. Therefore, the CAPRI method was used in this study. The CAPRI method was put forward by Piscini et al. (2017) to analyse the Amatrice-Norcia seismic sequence (Central Italy), which began on August 24, 2016 with an M6 major shock [26]. The same approach was used to evaluate the effectiveness of the method for detecting anomalous signals for impending volcanic eruptions [39].

The method is based on the comparison of about forty-year historical series in the same seasonal period of land/atmospheric parameters and the year in which the seismic sequence occurred. The approach lies in the way in which the time series are processed, in particular by identifying anomalous values with respect to the normal historical background and correctly eliminating the possible effect of global warming. The description of CAPRI method is as below:

where $T_{ERA}(d)$ is the mean value of SKT on the d day (13 Feb–12 May) in the year (1980–2007); λ, φ are the latitude and longitude, respectively; triangular brackets represent the spatial average; m(d) is the fit slope; $\Delta T(d)=\mathrm{m}\left(\mathrm{d}\right)\times {(\mathrm{y}-\mathrm{y}}_0)$ is then used to remove the long-term variation of the variable analysed for the same day (d); y is the considered year and ${\mathrm{y}}_0$ is the first year of time series and is used as a “reference”. $T_{\mathrm{h}}(d)$ is the average temperature of the particular day in the years prior to the EQ; $N_{\mathrm{y}}$ is the total number of analysed time series; in the process of calculating formula (3), the standard deviation σ is also calculated. The variable behaviour of the year in analysis ($\tilde{\mathrm{y}}$) is then compared with the historical series.

MERRA-2 atmospheric data were analyzed using an adapted version of CAPRI algorithm, where neither land-sea masking nor global warming removing procedure was applied.

4. Result

4.1. Case Study: Wenchuan EQ

4.1.1. SKT

Figure 2 shows the statistical comparison between the 2008 and the historical 1980-2007 time series for SKT: the time interval ranges from 13 February to 12 May, the first day and the 90st day indicated by the x-coordinate, respectively. The temperature is more than 2σ higher than the historical for three days (red circle in Figure 2 highlights this 3-day persistent anomaly): the persistence of this feature for three consecutive days (around 40 days before the mainshock) is noteworthy. Before and after these three days, there were also three individual points whose SKT values were higher than 2σ.

The spatial distribution of SKT anomalies in these days is shown in Figure 3. On the 28th day (12, March), SKT high value area first appeared near the left of the epicenter. On the 55th day, SKT high values were not around the epicenter. Then, for three days, from the 57th to 59th day, we observe that the SKT high value area was always located near the epicenter. On the 65th day (18, April), SKT high values located at the left of the epicenter of Wenchuan EQ. Therefore, we considered SKT anomaly phenomena of the three days (from the 57th to 59th) might be related to the Wenchuan EQ.

4.1.2. SLHF

Figure 4 describes the statistical comparison between the 2008 SLHF and the historical 1980–2007 time series: the x axis ranges from 13 February to 12 May. The SLHF is almost greater than 1.5σ from the historical mean for two days highlighted in Figure 4, almost simultaneously with SKT anomaly.

For the 27th and 54th days, the spatial distributions of SLHF parameter were extracted and shown in Figure 5, where we can notice that SLHF high-value areas of the 27th day were along the LMSF zone, and that on the 54th day the SLHF high-value area was located in the northeast of the epicenter.

4.1.3. TCWV

The statistical analysis for TCWV is shown in Figure 6, where the comparison between the 2008 TCWV and the historical 1980-2007 time series is shown: the x-axis ranges from 13 February to 12 May. The TCWV values exceed 2σ the historical mean for two single days.

The TCWV anomalous spatial location was extracted and the results were shown in Figure 7. High TCWV values could be observed near the LMSF zone and the east of the fault. The water vapor in the atmosphere comes from the underlying surface, including the evaporation of the water surface, the surface of wet objects, and the foliage of plants. The TCWV anomalous spatial location might be affected by geographical location and topographical factors.

4.1.4. AOD

We report the statistical analysis for AOD in Figure 8, where we show the comparison between the 2008 AOD and the historical 1980-2007 time series: the range starts from 13 February to 12. The AOD values are greater than 2σ than the historical mean for several days (the red circle in Figure 8 highlights AOD anomaly).

The AOD spatial distribution on days the 4–6th, 12th, 21–22th, 57–60th and 80th was shown in Figure 9. The results showed that the high AOD values were mainly distributed along the LMSF zone on the 4–6th days, and the main AOD distribution on the 21–22th days was first along the LMSF zone and then the northeast direction of the epicenter. The single-day high value distribution may be contrary to the continuity of seismic activity. The continuous high values of AOD on the 57th to 60th days are mainly distributed in the epicenter of Wenchuan EQ and the LMSF zone. On the 80th day (10 days before the EQ), AOD abnormal high values appeared around the epicenter, and the high AOD value was very obvious. AOD persistently high values were mainly distributed near the epicenter and along with the LMSF zone, which might be related to the tectonic activity of the LMSF zone.

4.2. Case Study: Lushan EQ

4.2.1. SKT

In Figure 10, the comparison between the 2013 SKT and the historical 1980-2012 time series is shown: the x axis ranges from 21 January to 20 April. The results showed that SKT on the 56th day (March 17) before the Lushan EQ was higher 2σ than the historical mean. Although there was only one SKT value higher 2σ than the historical mean, the SKT value on the previous day (March 16) was also high, almost approaching the 2σ (a difference of nearly 0.1 K) and almost have the trend of abnormal SKT increase lasting for two days.

Further, the high-value SKT spatial distribution on the 56th day (March 17) before the Lushan EQ was observed, as shown in Figure 11. SKT high value area was located near the epicenter of Lushan EQ and LMSF zone. Therefore, it was considered that SKT anomaly on March 17 was related to Lushan EQ.

4.2.2. SLHF

In Figure 12, we show the statistical comparison between the 2013 SLHF and the historical 1980-2012 time series: the x axis ranges from 21 January to 20 April. The results showed that SLHF values before the Lushan EQ were lower 2σ than the historical mean. Only on two single days, SLHF values were higher than 1.5σ than the historical mean.

The high-value SLHF spatial distribution on the 38th and 60th day was observed, as shown in Figure 13. It could be observed that the high-value area of SLHF is distributed near the LMSF zone with a strip.

4.2.3. TCWV

In Figure 14, we show the statistical comparison between the 2013 TCWV and the historical 1980-2012 time series: the x axis ranges from 21 January to 20 April. The results showed that TCWV values on the 17~18th, 28th and 55th days before the Lushan EQ were higher 2σ than the historical mean for three days.

Further, the high-value TCWV spatial distribution during four days before the Lushan EQ was observed, as shown in Figure 15. TCWV high value areas were located near the epicenter and the northeast of Lushan EQ on the 17th, 18th, 55th days, prior to Lushan EQ.

4.2.4. AOD

The statistical comparison between the 2013 AOD and the historical 1980-2012 time series is shown in Figure 16, where the x axis ranges from 21 January to 20 April. The results showed that AOD values on the 7th (January 27), 18th (February 15), 49~53th (March 9~13), 78~80th (April 7~9) days before the Lushan EQ were higher 2σ than the historical mean for these days.

Further, the high-value AOD spatial distributions before the Lushan EQ were observed, as shown in Figure 17 and Figure 18. On the 7th day, the AOD high value was obviously distributed in the northeast direction of the epicenter, and the AOD was high along the LMSF zone. On the 18th day, the AOD high value area was not obvious; on the 49-53th days, the AOD anomalous high value area was mainly located in the epicenter and Longmeshan fault. The time evolution of the fault zone was that it first strengthened and then weakened along the LMSF zone. The high value of AOD on the 52th day is mainly distributed near the epicenter and weakened on the 54th day. The AOD difference results on the 78th to 80th days are not obvious and the high value area mainly distributed southwest of the epicenter. The AOD results showed that the evolution process of AOD along the LMSF zone from the 49th to 53th days is obvious.

5. Discussion

The time sequence of the four atmospheric parameters for Wenchuan and Lushan EQs is summarized in Figure 19. Despite the wide time span, the anomaly times of four parameters are focused on two months before Wenchuan and Lushan EQs. In terms of the time sequence of the anomalous behaviour of the four parameters, AOD and TCWV responses are the earliest, followed by SKT and SLHF in both Wenchuan and Lushan EQs. First of all, both EQs were located along the LMSFsystem. The maps of the parameters SKT (Figure 11), AOD (Figure 17 and Figure 18), SLHF (Figure 13) clearly show the close spatial correlation between the anomalies and the location of the faults responsible for the Wenchuan and Lushan EQs. This spatial correlation is less evident for TCWV (Figure 7), whose anomaly is, however, close to the fault segment responsible for the Lushan event. This strict correspondence between the LMSF system and the parameters analyzed here indicates that the observed anomalies are related to some processes occurring near/on the faults. Seismological evidences of foreshocks before the Wenchuan and Lushan EQs are lacking [41,42]. This suggests that the anomalies related to the SKT, AOD, SLHF, and TCWV are not related to stress changes induced by tectonic processes, i.e., early fracturation. Alternatively, an aseismic mechanism of diffuse fracturation by stress changes within a ductile crust may be invoked; however, the Wenchuan and Lushan mainshocks and their aftershocks occurred within two high Vp and Vs zones of the 18-20 km deep Guanxian-Anxian Fault and Yingxiu-Beichuan Fault [34]. Therefore, a mechanism of ductile deformation with diffuse fracturing is not easy to be responsible for the observed anomalies. Tomographic images [43] and magnetotelluric data [36] show that two low velocity zones separating the focal volumes of the Wenchuan and Lushan mainshocks and aftershocks extend in depth from the middle-crust to the surface. These zones are interpreted as fractured, fluid-rich crustal volumes. As already reported, Zhao et al. (2012), evidence a low-resistivity crustal zone close to the LMSF system [37]. This zone, which is located at the base of the upper crust is interpreted as fluid-rich, mechanically weak layer. Lei and Zhao (2009) propose that these deep fluids may be responsible for the fault weakening of the LMSF systems [38]. Data on gas emissions by Cui et al., (2017; 2019) show huge anomalies in the CH₄ and CO₂ discharge from the LMSF system before the Lushan and Wenchuan EQs [39,44]. These data and the occurrence of peak values in the detected anomaly at the intersection area of faults provide a clear evidence of the relative larger permeability of the LMSF system with respect to the surrounding Tibetan plateau and Sichuan basin (Figure 1). According to Cui et al. (2017), the gas emissions along the faults before the two EQs may induce a variety of anomalies in the atmosphere. On the basis of our data and the above described crustal structure of the LMSF system and its degassing processes, we propose that the anomalies we detect in the SKT, AOD, SLHF, and TCWV parameters before the Lushan and Wenchuan events may be due to the anomalous discharge of CH₄ and CO₂ along the seismogenic faults. In this framework, the passive upraise of fluids from the middle/upper crust to the surface may modify the elastic properties of rocks (fault weakening), saturate the local fracture system and induce stress variations [7,29,45]. We exclude an active upraising of over-pressurized fluids within the faults because such fluids should have triggered EQs, which is not observed before Lushan and Wenchuan. Our conceptual model is consistent with the lithosphere and atmosphere coupling model (LAIC). Changes in crustal stress due to fluid movements as well as the upraising of hot fluids before Wenchuan and Lushan EQs may induce thermal anomalies. AOD anomaly has been reported in Wenchuan and Lushan EQs by other methods [16,30] and the AOD anomaly response time we found here was similar. SKT variations were most likely caused by both geological processes and action of fluid related crustal stress variations in the crust, which resulted in the Wenchuan and Lushan EQs [38]. Our data show that the SLHF anomaly almost simultaneously occur with the TCWV anomaly. This result agrees with the model of Dey et al. (2003, 2004), which explains the increase in evaporation of the surface water because of the increased SLHF [12,46]. Therefore, when a large amount of fluids possibly associated to stress variations is released to the atmosphere, water vapor formation is promoted. These processes may cause the anomalies of water vapor distributed along the seismic fault zones. The SLHF anomaly detected before the Wenchuan event is probably related to a variation in surface temperature promoted by the release of deep CO₂ [39]. The involvement of deep fluids we propose for the studied compressional-type Wenchuan and Lushan events has been also recognised in strike-slip [47,48] and extensional [49,50,51] EQs and seismic sequences.

From the perspective of spatial location, the high values of SKT parameter for Wenchuan and Lushan are located in the epicentral areas. The AOD anomaly values concentrate around the LMSF system. The SLHF abnormal area of Wenchuan event is also located in the LMSF system and its northeast region. The SLHF abnormal area associated to the Lushan event describe a strip along the LMSF zone. TCWV associated to Wenchuan is also located along a strip around the LMSF zone. TCWV anomaly of Lushan is on the northeast side of the LMSF zone. Previous studies have suggested that the Wenchuan event occurred in the central LMSF zone, and that the Lushan event was located at the eastern side of the fore fault of the southwest segment of LMSF zone [29,32,52]. It was also reported that fluid-bearing ductile flow was from the lower crustal materials of Tibet being pushed into a weakened segment of the LMSF zone. Therefore, the four-parameter anomaly areas we report might be related to the processes preparing the fault rupture. Xie et al. (2015) showed that the relative wavelet power spectrum (RWPS) anomalies of brightness temperature data of the Lushan event firstly appeared at the middle part of LMSF zone [53]. This is in accordance with our findings on the location of the SKT, SLHF and TCWV anomaly zones. From the four parametric anomaly information areas, the Lushan seismic anomaly response area includes the LMSF zone. It was reported by Li et al. (2018) that the 2013 Ms 7.0 Lushan earthquake should be defined as a “delayed” aftershock of the 2008 Ms 8.0 Wenchuan earthquake by the quantity analysis on the seismic moment accumulated and released along the southern LMSF zone [54]. Furthermore, the calculated static Coulomb stress changes also indicated that the 2013 Lushan earthquake might have been statically triggered by the 2008 Wenchuan earthquake [25]. Therefore, the anomalies around the LMSF zone related to the Lushan EQ could be partially attributed to the Wenchuan EQ because of the seismotectonic changes and fractures in the crust produced by the Wenchuan EQ.

6. Conclusions

The all anomalous responses of four climatological parameters occurred within the two months before Wenchuan and Lushan EQs. In terms of the time response ordering of the four analyzed parameters, AOD and TCWV responses are the earliest, followed by SKT and SLHF in both Wenchuan and Lushan EQs. Both EQs were located in the LMSF system, and the similarity of time sequence of four parameters may be related to the processes occurring on the active LMSF segments. The anomalies of water vapor distributed along the seismic fault zones were caused by the upraising and release to the atmosphere of deep, crustal fluids during the preparation phase of the EQs [46]. AOD anomaly may be related to SLHF and TCWV and SKT. From the perspective of spatial location, the larger values of the anomalies of the four parameters are located in the Wenchuan and Lushan epicentral areas and along the LMSF zone. This result is consistent with both fluid migration processes and crustal stress changes. The anomalies related to the Lushan EQ could be partially attributed to the Wenchuan EQ because of the new fractures and increased permeability in the crust produced by the Wenchuan event. The results of this study indicate that the multi-parametric approach is beneficial to capture anomalous phenomena before strong EQs and highlight the role of fluids in the preparing phase of large EQ.