Citation: | Wang, H., Song, D. M., Shan, X. J., and Wang, B. (2024). A method for extracting the preseismic gravity anomalies over the Tibetan Plateau based on the maximum shear strain using GRACE data. Earth Planet. Phys., 8(4), 589–608. DOI: 10.26464/epp2024023 |
The occurrence of earthquakes is closely related to the crustal geotectonic movement and the migration of mass, which consequently cause changes in gravity. The Gravity Recovery And Climate Experiment (GRACE) satellite data can be used to detect gravity changes associated with large earthquakes. However, previous GRACE satellite-based seismic gravity-change studies have focused more on coseismic gravity changes than on preseismic gravity changes. Moreover, the noise of the north–south stripe in GRACE data is difficult to eliminate, thereby resulting in the loss of some gravity information related to tectonic activities. To explore the preseismic gravity anomalies in a more refined way, we first propose a method of characterizing gravity variation based on the maximum shear strain of gravity, inspired by the concept of crustal strain. The offset index method is then adopted to describe the gravity anomalies, and the spatial and temporal characteristics of gravity anomalies before earthquakes are analyzed at the scales of the fault zone and plate, respectively. In this work, experiments are carried out on the Tibetan Plateau and its surrounding areas, and the following findings are obtained: First, from the observation scale of the fault zone, we detect the occurrence of large-area gravity anomalies near the epicenter, oftentimes about half a year before an earthquake, and these anomalies were distributed along the fault zone. Second, from the observation scale of the plate, we find that when an earthquake occurred on the Tibetan Plateau, a large number of gravity anomalies also occurred at the boundary of the Tibetan Plateau and the Indian Plate. Moreover, the aforementioned experiments confirm that the proposed method can successfully capture the preseismic gravity anomalies of large earthquakes with a magnitude of less than 8, which suggests a new idea for the application of gravity satellite data to earthquake research.
An earthquake can be an extremely destructive geological disaster, causing great loss of lives and property and the disruption of the national economic infrastructure (Allen and Melgar, 2019; Fan XM et al., 2019). Therefore, seismic monitoring is of great practical significance for seismic fortification (Tan CX et al., 2019). For this reason, several seismic monitoring stations have been established in earthquake-prone countries or regions around the world (Saç et al., 2011; D’Alessandro et al., 2018). However, traditional means of seismic monitoring, such as geodetic networks, can record seismic precursors only within a limited area, and the seismogenic process is difficult to replicate on a macro level (Zieger and Ritter, 2018; Martinelli, 2020).
In the late 20th century, the rapid development of remote sensing technology enabled the realization of large-area observation and real-time information acquisition, which effectively made up for the shortcomings of traditional means of seismic monitoring, such as fixed-point observation (Marchetti et al., 2020; Mahmood, 2021; Zhao XW et al., 2021). For example, satellite radar interferometry could be used to obtain the postseismic deformation field covering the entire seismic fault zone (Massonnet et al., 1994; Brozzetti et al., 2020). Qu W et al. (2018) studied present-day crustal deformation characteristics of the southeastern Tibetan Plateau and the surrounding areas based on global positioning system (GPS) data for the periods from 1999 to 2007 and 2009 to 2011. They concluded that regions with frequent earthquakes often correspond to the transition zone between regions with high shear strain and a plane strain gradient. Although the interferometric synthetic aperture radar (InSAR) technology and GPS measurements mentioned can better monitor the crustal movement on the surface, it is difficult to detect the material migration and mass change in the deep crust (Panet et al., 2010).
The gravity field reflects the internal structure and composition of the Earth and is regarded as a historical representation of the dynamic processes of the solid Earth (Sun HP et al., 2017; Sansò and Barzaghi, 2020; Duan HR et al., 2022). Before and after the 1964 Alaska earthquake, the 1965–1967 Songdai earthquake swarm in Japan, the 1974 Izu earthquake in Japan, the 1976 Gazley earthquake swarm in the Soviet Union, and the 2013 Ms 7.0 Lushan earthquake and gravity drop after the mainshock (Barnes, 1966; Wang LH et al., 2023 Zhu YQ et al., 2023), the time-varying information on gravity before earthquakes was observed with time series of absolute and relative gravity data. The gravity data obtained by the Gravity Recovery And Climate Experiment (GRACE), which is representative of gravity satellite technology, contains information on the crustal density distribution and can be used to study the interaction mechanism between and within plates before earthquakes (Wahr et al., 1998; Peidou and Pagiatakis, 2019; Zotin et al., 2019; Gido et al., 2020; Velicogna et al., 2020; Liu C and Sun WK, 2023). Sun W and Okubo (2004) studied whether the GRACE satellites could detect coseismic gravity changes based on dislocation theory and concluded that GRACE is capable of detecting coseismic deformation signals from earthquakes with a magnitude greater than 7.5.
Since then, gravity signals caused by several large earthquakes, such as the 2004 Sumatra–Andaman earthquake, the 2005 MW 8.6 Indonesia earthquake, the 2010 MW 8.8 Chile earthquake, and the 2011 MW 9.0 Japan earthquake, have successfully been extracted by using GRACE data (Han SC et al., 2006; Panet et al., 2018; Chao BF and Liau, 2019; Mikhailov et al., 2020). However, these studies have mainly focused on the dynamics of the coseismic gravity field rather than on preseismic gravity anomalies, which are of great significance for the early warning of earthquakes (Sasmal et al., 2021). For example, Chen SM et al. (2016) found a significant gravity increase, up to 20 μGal/year, in southern Tibet between 2010–2011 and 2013. They suggested that these gravity changes may have been related to strain accumulation and possibly mass migration in a broad source region of the 2015 Nepal earthquake. After the Wenchuan MS 8.0 earthquake in 2008, the integration of the national network and the regional network accelerated, forming the whole gravity observation network on the Chinese continent. This network has made relatively successful medium-term predictions for a series of earthquakes with MS 6.0 or above (such as the MS 7.0 Lushan earthquake, the MS 6.4 Menyuan earthquake, and the MS 7.0 Jiuzhaigou earthquake) in recent years and has played an important role in the study of earthquake mechanisms and earthquake prediction levels in China (Zhu YQ et al., 2023). Moreover, when obtaining gravity anomalies, the commonly adopted difference disposal of the gravity field with the gravity field of adjacent months or the average gravity field of many years cannot effectively remove the inherent north–south stripe noise in GRACE data (Wang XL et al., 2019; Peidou and Pagiatakis, 2020). On the contrary, it is more likely to cause the annihilation of medium- and high-order information in the GRACE gravity field model, which results in the loss of some gravity information related to tectonic activity.
To explore preseismic gravity anomalies in a more refined way, in this study we propose a method of characterizing gravity variation based on the maximum shear strain of gravity, inspired by the concept of crustal strain (Dermanis and Livireratos, 1983). The offset index method (Tramutoli, 1998; Song DM et al., 2018) is then adopted to describe the gravity anomalies existing in the time series of the maximum shear strain. The experimental results on the Tibetan Plateau and surrounding areas with frequent strong earthquakes confirm that the proposed method can successfully extract the preseismic gravity anomalies of large earthquakes with a magnitude of less than 8.
The study area selected in this research is the Tibetan Plateau and its surrounding areas. Recognized as the highest, largest, and most complex plateau in the world, the Tibetan Plateau is located in the middle of the Asian continent and extends over the latitude–longitude domain of 26°00′−39°47′N, 73°19′−104°47′E. Since the Cenozoic era, collision between the Indian Plate and the Eurasian Plate has formed a large-scale crustal uplift of the Tibetan Plateau, resulting in a crust approximately 80 km in thickness (Cui et al., 2003; Figure 1a). Because of warming of the crust, the internal medium is highly viscoplastic, and the Indian Plate continues to push the plateau northward. This causes the lower crust to constantly drive the movement and deformation of the brittle layer in the upper crust, thus further forming the complex seismogenic environment of the Tibetan Plateau (Royden et al., 2008), which is the main area where strong land earthquakes occur globally.
Since the 21st century, the tectonic movement of the Tibetan Plateau has typically been active, and strong earthquakes have occurred frequently (Molnar and Lyon-Caent, 1989). From 2004 to 2016, 32 strong earthquakes of magnitude 6 or above were counted in this region, including two large earthquakes of magnitudes greater than 7.5, namely the MW 7.9 Wenchuan earthquake and the MW 7.8 Nepal earthquake, which occurred in the Longmenshan Fault Zone and the southern Himalayan Orogenic Belt, respectively (U.S. Geological Survey, http://earthquake.usgs.gov/earthquakes/search/). Specifically, the Longmenshan Fault Zone (F1 in Figure 1b), located in the eastern part of the Tibetan Plateau, is a NE–SW-extending structural belt with a length of more than 500 km and a width of 30−50 km (Dirks et al., 1994). In addition, the Main Himalayan Thrust Belt (F2 in Figure 1b) is a southward-protruding E–W arc tectonic belt, approximately 2500 km long and 300−500 km wide, which belongs to the southern boundary of the plateau and is situated at the junction of the Indian Plate and the Eurasian Plate (Burg et al., 1997). It should be noted that within 3 months after the Wenchuan and Nepal earthquakes, several earthquakes of magnitude 6 or above occurred in the seismogenic fault zone. Detailed information on these earthquakes is given in Table 1.
Seismogenic fault | Data | Epicenter location | MW | Depth (km) | |
Longitude (°E) | Latitude (°N) | ||||
Longmenshan | 2008 May 12 | 103.322 | 31.002 | 7.9 | 19.0 |
2008 May 12 | 103.618 | 31.214 | 6.1 | 10.0 | |
2008 May 25 | 105.423 | 32.56 | 6.1 | 18.0 | |
2008 August 12 | 105.494 | 32.756 | 6.0 | 6.0 | |
Himalayan Orogenic Belt | 2015 April 25 | 84.731 | 28.231 | 7.8 | 8.2 |
2015 April 25 | 85.540 | 27.629 | 6.1 | 10.0 | |
2015 April 25 | 84.822 | 28.224 | 6.6 | 10.0 | |
2015 April 26 | 86.017 | 27.771 | 6.7 | 22.9 | |
2015 May 12 | 86.066 | 27.809 | 7.3 | 15.0 | |
2015 May 12 | 86.162 | 27.625 | 6.3 | 15.0 |
The GRACE satellites, jointly designed by the Deutsches Zentrum für Luft- und Raumfahrt (DLR) and the National Aeronautics and Space Administration (NASA), was successfully launched in March 2002 (Ince et al., 2019). The GRACE satellite system consists of two low-orbit satellites that operate in satellite tracking mode, which includes a low-Low (SST-LL) mode and a high-Low (SST-HL) mode. We use GRACE Release 6 (RL06) of the Level-2 product in this work. This product was released by the Center for Space Research (CSR) at the University of Texas at Austin (http://www2.csr.utexas.edu/grace/RL06.html). It contains 163 months of gravitational field data from April 2002 to July 2017, with a maximum order of 96, a spatial resolution of 300 km, and a temporal resolution of one month. Note that this product excludes the tidal and nontidal variations in atmospheric and ocean mass (Swenson et al., 2008). Because GRACE is a polar-orbiting satellite with a large orbital inclination, the accuracy of the degree 2 coefficients of the GRACE gravity field model is low. In this work, the degree 2 coefficients are thus replaced by the measurement values of satellite laser ranging (SLR) provided by the CSR (Cheng MK and Tapley, 2004).
The Global Land Data Assimilation System (GLDAS) is a joint project of NASA, the National Centers for Environmental Prediction (NCEP), and the National Oceanic and Atmospheric Administration (NOAA; Rodell et al., 2004) using data assimilation technology to integrate satellite and surface observations into a unified model. The various land surface field information is organized in the form of grid data, such as global precipitation (rainfall and snowfall), evapotranspiration (soil water evaporation and vegetation transpiration), surface runoff, subsurface runoff, soil moisture, surface temperature, and surface heat flow. The GLDAS-2.1 model data (https://ldas.gsfc.nasa.gov/gldas/) include soil moisture 0–2 m underground and snow water, with a temporal resolution of one month and a spatial resolution of 1°.
To explore information on the tectonic activity before earthquakes, here we propose a method for extracting gravity anomalies on the basis of the maximum shear strain of gravity. Because the maximum shear strain is calculated by the second-order gradient of the GRACE disturbance potential, this method can suppress the stripe noise better than can difference disposal, thus effectively improving the sensitivity of gravity anomaly detection. A detailed description of the process is given in Figure 2.
This method consists of three specific steps:
(1) Removal of the effects of seasonal changes in hydrologic information. To extract gravity signals related to tectonic activity, the periodic changes in the gravity field caused by seasonal variations in hydrology should be removed first so as to obtain the increment of disturbance potential without interference of the hydrology. See Section 3.1 for details.
(2) Calculation of the maximum shear strain of the crust before earthquakes. After removing the influence of hydrology, gravity information contains much purer information on crustal movement and material migration, but it is still necessary to pick out some physical parameters to achieve the quantitative characterization of pre-earthquake tectonic activities. Inspired by the concept of crustal strain, we introduce the maximum shear strain of gravity in this study. In other words, on the basis of Step 1, the gravity strain tensor is obtained by further calculating the second-order gradient of the increment of disturbance potential, and then the maximum shear strain of gravity is ultimately generated to characterize the pre-earthquake tectonic activities. Simultaneously, note that the stripe noise is effectively suppressed once the second-order gradient of the disturbance potential is used. The specific procedures are described in Section 3.2.
(3) Spatiotemporal analysis of anomaly information before the earthquake. To further characterize the anomalies in pre-earthquake tectonic activities and explore the spatiotemporal evolution of pre-earthquake anomalies, the offset index K is calculated on the basis of the time series of maximum shear strain. Furthermore, the spatial distribution characteristics of the K value are analyzed. Section 3.3 is devoted to a detailed description of the implementation process.
The gravity changes can be recovered based on GRACE satellite data in terms of spherical harmonic coefficients, which are actually the solution to the partial differential equation of the disturbance potential, satisfying the following Laplace equation (Sun W and Okubo, 2004; Eshagh and Abdollahzadeh, 2012; Tang H and Sun W, 2021):
T(r,θ,φ)=GMR∞∑n=2(Rr)n+1n∑m=−n(Cnmcosmφ+Snmsinmφ)ˉPnm(cosθ), | (1) |
where GM is the geocentric gravitational constant, (r, θ, φ) represents the spherical coordinates, R is the mean radius of the Earth,
In addition, the potential coefficient of the gravitational field can be used to calculate the density change of the Earth’s surface so as to represent the spatial distribution change of the Earth’s surface mass, which can be expressed as
Δσ(θ,φ)=2RρEπ3∞∑n=22n+11+knn∑m=−n(Cnmcosmφ+Snmsinmφ)ˉPnm(cosθ), | (2) |
where,
However, the changes in surface mass calculated by GRACE data are caused by the combination of tectonic activities and seasonal changes in hydrology (Rao W and Sun WK, 2022). Furthermore, the gravity variations caused by seasonal changes in terrestrial water storage (e.g., abundant rainfall in the Sichuan Basin in summer; Tang L et al., 2020) are relatively significant. In particular, the influence of hydrological signals in GRACE data can be eliminated by utilizing the soil water and snow water data provided by the GLDAS hydrological model. For this purpose, the spherical harmonic coefficients of GLDAS data are first expanded and cut off to degree 96 (consistent with GRACE observations) to obtain the hydrological spherical harmonic coefficients Cgldas and Sgldas. After that, by the difference between Cgldas, Sgldas and the spherical harmonic coefficients Cgrace, Sgrace of GRACE data, respectively, the corrected spherical harmonic coefficients
{ΔCti=Ccorrectti+1−Ccorrectti,ΔSti=Scorrectti+1−Scorrectti, | (3) |
where
Let
Gravity changes refer to the variations in gravity at some point on the Earth’s crust at two distinct time states (Dermanis and Livireratos, 1983; Fatolazadeh et al., 2020). In this study, considering that the temporal resolution of GRACE data is one month, we hereby use adjacent months as observation windows of gravity field deformation, where
{{\boldsymbol{E}}_i} = \frac{1}{2}\left[ {{{\left( {\frac{{\partial {\text{γ}} }}{{\partial {\text{Γ}} }}} \right)}^{\mathrm{T}}}\left( {\frac{{\partial {\text{γ}} }}{{\partial {\text{Γ}} }}} \right) - {\boldsymbol{I}}} \right] = \frac{1}{2}\left( {{\boldsymbol{B}}_i^{ - 1}{\boldsymbol{B}}_{i + 1}^2{\boldsymbol{B}}_i^{ - 1} - {\boldsymbol{I}}} \right) , | (4) |
where
{\text{Γ}}=\dfrac{\partial T_{i}(\boldsymbol {R})}{\partial \boldsymbol {R}} , {\text{γ}}=\dfrac{\partial T_{i+1}(\boldsymbol {R})}{\partial \boldsymbol {R}}, |
and
{{\boldsymbol{B}}_i} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{{\partial ^2}{T_i}}}{{\partial {X^2}}}}&{\dfrac{{{\partial ^2}{T_i}}}{{\partial X\partial Y}}}&{\dfrac{{{\partial ^2}{T_i}}}{{\partial X\partial Z}}} \\ {\dfrac{{{\partial ^2}{T_i}}}{{\partial Y\partial X}}}&{\dfrac{{{\partial ^2}{T_i}}}{{\partial {Y^2}}}}&{\dfrac{{{\partial ^2}{T_i}}}{{\partial Y\partial Z}}} \\ {\dfrac{{{\partial ^2}{T_i}}}{{\partial Z\partial X}}}&{\dfrac{{{\partial ^2}{T_i}}}{{\partial Z\partial Y}}}&{\dfrac{{{\partial ^2}{T_i}}}{{\partial {Z^2}}}} \end{array}} \right] . | (5) |
Actually, the Lagrange gravity strain tensor calculated by Equation (5) is a symmetric matrix. In general, ε11, ε12, and ε13 represent the components of elongation in the three coordinate directions:
{\boldsymbol{E }}= \left[ {\begin{array}{*{20}{c}} {{\varepsilon _{{\text{11}}}}}&{{\varepsilon _{{\text{12}}}}}&{{\varepsilon _{{\text{13}}}}} \\ {{\varepsilon _{{\text{21}}}}}&{{\varepsilon _{{\text{22}}}}}&{{\varepsilon _{{\text{23}}}}} \\ {{\varepsilon _{{\text{31}}}}}&{{\varepsilon _{{\text{32}}}}}&{{\varepsilon _{{\text{33}}}}} \end{array}} \right] . | (6) |
The quantities λmax, λmin, which are the maximum and minimum eigenvalues of the gravity strain tensor E above, are called the (Lagrangian) principal strains. Indeed, the principal strains λmax, λmin are the eigenvalues of E and can be computed as the roots of the characteristic equation:
\mathrm{det}\left(\boldsymbol{E}-\lambda\boldsymbol{I}\right)=\lambda^2-\mathrm{tr}\left(\boldsymbol{E}\right)\lambda+\mathrm{det}\left(\boldsymbol{E}\right)=0, | (7) |
which are
{{\lambda }}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{1}{2}\left\{{\mathrm{tr}}\left({\boldsymbol{E}}\right)+{\left[{{\mathrm{tr}}\left({\boldsymbol{E}}\right)}^{2}-4\mathrm{d}\mathrm{e}\mathrm{t}\left({\boldsymbol{E}}\right)\right]}^{1/2}\right\} , |
{{\lambda }}_{\mathrm{m}\mathrm{i}\mathrm{n}}=\frac{1}{2}\left\{{\mathrm{tr}}\left({\boldsymbol{E}}\right)-{\left[{{\mathrm{tr}}\left({\boldsymbol{E}}\right)}^{2}-4\mathrm{d}\mathrm{e}\mathrm{t}\left({\boldsymbol{E}}\right)\right]}^{1/2}\right\} . | (8) |
The maximum shear strain (MSH) is defined as the difference between the maximum and minimum eigenvalues of the gravity strain tensor E above, representing the anisotropic part of deformation in the infinitesimal vicinity of the point of interest:
\mathrm{M}\mathrm{S}\mathrm{H}={\lambda}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\lambda}_{\mathrm{m}\mathrm{i}\mathrm{n}} . | (9) |
The physical interpretation of MSH is the shear strain across the direction of its maximum value (always positive) under the action of stress. For this reason, we use the maximum shear strain of gravity in this study as a characterization to describe the strength of tectonic activity.
To better understand the seismogenic process of the fault zone by further extracting the pre-earthquake anomalous changes, data on the maximum shear strain time series are analyzed in this study by means of the offset index K (Tramutoli, 1998; Song DM et al., 2018). Additionally, the threshold of the offset index K value needs to be preset for anomaly detection of different time series at each pixel point. The calculation formula of offset index K is given as follows:
K(x,y,t)=\frac{\mathrm{M}\mathrm{S}\mathrm{H}(x,y,t)-{\mu }_{\mathrm{M}\mathrm{S}\mathrm{H}}(x,y)}{{\sigma }_{\mathrm{M}\mathrm{S}\mathrm{H}}(x,y)} , | (10) |
where
Figure 3 shows the process of calculating the offset index K based on the maximum shear strain of gravity. According to the Pauta criterion of statistics, a sample is judged as an outlier if it falls outside the range of three standard deviations from the mean of the sample set (Li LM et al., 2016; Seheult et al., 1989). This is the basis for selecting the threshold value of offset index K as 3 in this study. Note that the calculation of the offset index value is aimed at changes in each pixel in its time series and has nothing to do with other pixels in the spatial domain. In addition, from detecting the time series anomalies in each pixel, the distribution characteristics of the K value of each pixel in the spatial domain are also observed so as to obtain an understanding of the spatiotemporal evolution of anomalies.
In addition, the scope of seismogenic zones corresponding to different magnitudes is different. Dobrovolsky et al. (1979) gave the formula for calculating the seismogenic area, based on a large number of experiments and observation results, as shown in Equation (11):
R = {10^{0.43M}} , | (11) |
where R is the radius of the seismogenic area and M is the earthquake magnitude. According to the equation above, the larger is the earthquake magnitude, the larger will be the corresponding seismogenic area.
In this work, the Tibetan Plateau and its surroundings are taken as the study area. The gravity strain tensor is calculated according to the second-order gradient of the increment of the disturbance potential. In addition, the time series of the maximum shear strain of gravity is derived for each pixel 2 years before the earthquake. The gravity anomalies in the time series can then be depicted by the offset index K. Finally, the temporal and spatial variation characteristics of gravity anomalies before earthquakes are analyzed on the observation scale of the fault zone and the Tibetan Plateau, respectively, to explore the possible pre-earthquake tectonic activities in the Longmenshan Fault Zone, the Himalayan Orogenic Belt, and the Tibetan Plateau.
The occurrence of earthquakes is closely related to the tectonic activities in fault zones. To further explore the seismogenic mechanism, in this study we investigate the characteristics of pre-earthquake tectonic activity of the Longmenshan Fault and the southern Himalayan Orogenic Belt, which are the seismogenic fault zones responsible for the Wenchuan earthquake and the Nepal earthquake, respectively. In addition, to verify that the gravity anomaly is a unique phenomenon before earthquakes, comparison experiments are carried out between nonearthquake years and seismogenic years. In particular, the offset index values for the nonearthquake years are calculated by using GRACE satellite data from 2005−2006, 2007−2008, and 2009−2010. The detailed calculation process is as follows. First, on the basis of the maximum shear strain of gravity, the offset index K values are calculated for each 0.5° × 0.5° grid within the fault zone area in the three periods mentioned above. The average of these three K values, in notation of K0, is then used to represent the corresponding state in the nonearthquake years. It should be noted that the nonearthquake years mentioned in this study do not specifically refer to years without earthquakes but to the time periods for comparison with the earthquake periods.
Because the GRACE observations also contain information on hydrological changes in the Earth’s mass, it is necessary to remove the hydrological anomalies and obtain information on the gravity anomalies caused only by tectonic activities before using the GRACE observations to study the gravity anomalies of earthquakes. To further illustrate the interference of hydrological information with the seismic signal, gravity change maps have been drawn showing the hydrological influence.
Figures 4a and 4b show the hydrological gravity distribution in the three fault zones studied in this work, with hydrological gravity values ranging from −1.21 to 1.30 and 0.70 to 6.27 μGal, respectively. The maximum hydrological gravity values within each region during the earthquake month, with 6.27 μGal values in the southern Himalayan Orogenic Belt region. However, Wang L et al. (2012) and Matsuo and Heki (2011) pointed out that the maximum gravity change caused by the largest earthquakes (M 9.2 Sumatra earthquake in 2004 and the M 9.0 Tohoku-oki earthquake in 2011) in the 21st century was approximately 10 μGal. Thus, this inconsistency implies that the interference of hydrological signals cannot be ignored in seismic research.
It can be further observed from the elliptical shapes in Figures 4c and 4d that, by comparison with the original data, after removing the hydrological signals, GRACE data were more stable during the nonseismic period, whereas the gravity change amplitude increased during the earthquake. For example, Figure 4c shows an obvious step-shaped decline in gravity value after removing the hydrological interference during the preseismic period. These bits of evidence indicate that the gravity anomaly information related to the tectonic activity of the fault zone is well highlighted after filtering out the hydrological factors from GRACE data, which further illustrates the importance of removing the hydrological interference.
Considering that the MW 7.9 Wenchuan earthquake and three aftershocks, namely, the MW 6.1 Wenchuan earthquake, the MW 6.1 Qingchuan earthquake, and the MW 6.0 Qingchuan earthquake, all occurred in the Longmenshan Fault Zone within 3 months, they are regarded as one earthquake sequence. Detailed information on these four earthquakes is shown in Table 1. Figure 5 shows the variation curves of total gravity offset index values (K > 0) and total gravity anomaly values (K > 3) in the Longmenshan Fault Zone 2 years before the earthquake. We find that 3 months before the Wenchuan earthquake, that is, half a year before the earthquake sequence, the gravity offset index values and gravity anomaly values in the fault zone both underwent a rapid increase simultaneously, and they reached their maximum value in 2 years. Therefore, it can be speculated that the tectonic activity in the fault zone was particularly active from February to April in 2008. To obtain the evolution law of gravity anomalies more clearly and intuitively before the earthquake, we obtain the spatiotemporal variations of the offset index K of the earthquake sequence that occurred in the Longmenshan Fault Zone half a year before the earthquake (Figure 6). The maximum K value is as high as 4.6, marked by the blue box in Figure 6b. Further observation shows that the number of high offset index values in Figures 6b and 6d is greater than those in other periods and is the most obvious in Figure 6b. Therefore, this result confirms the phenomenon that the tectonic activity in the fault zone was particularly active during the 3 months before the earthquake. For a better visual effect of the spatiotemporal changes in gravity anomalies before earthquakes, we specifically display the spatiotemporal distribution of K values greater than 3 in Figure 7.
Figure 7b shows that a large area of anomalies appears in the middle and north sections of the Longmenshan Fault Zone at 1 and 3 months before the MW 7.9 Wenchuan earthquake and its two aftershocks (the MW 6.1 Wenchuan earthquake and the MW 6.1 Qingchuan earthquake), that is, 4−6 months before the MW 6.0 Qingchuan earthquake, which is three times the area of the previous period (Figure 7a). This result shows that gravity anomalies had occurred in the middle and northern sections of the Longmenshan Fault Zone at this time. Combined with the fact that the GRACE gravity variation values calculated in this study were positive during this period, it can be inferred that stress accumulation had occurred on the fault, with a large amount of strain energy where the instability had reached a high level. Before this, only a few anomalies occurred in the southern part of the Longmenshan Fault Zone, as shown in Figure 7a. Note that the epicenters of the earthquake sequence are all located in the area of the anomaly mentioned above, among which the epicenters of the MW 6.1 and MW 6.0 Qingchuan earthquakes are located within 50 km of the maximum K value.
Figure 7d shows the spatial distribution of gravity anomalies one month after the MW 7.9 Wenchuan earthquake. We find that obvious gravity anomalies occurred again in the Longmenshan Fault Zone. In particular, high anomalies were distributed along the strike of the fault zone and mostly gathered around the north of the fault zone. According to the spatial distribution of gravity anomalies in Figures 7d and 7e, it can be inferred that the main tectonic activity zone of Longmenshan Fault Zone had moved northward at this time. Coincidentally, the MW 6.0 Qingchuan earthquake verifies this conjecture well because it occurred at the northeast end of the Longmenshan Fault Zone 2 months later. It also indicates that compared with a single earthquake, the spatiotemporal analysis of the earthquake sequence is able to capture characteristics of the tectonic activity in the fault zone more comprehensively. Figure 8 shows that during the same period of nonearthquake years, the maximum K value on the Longmenshan Fault Zone was only 2.8, implying that there was no gravity anomaly.
The Nepal earthquake sequence consists of six strong earthquakes that occurred less than a month apart, as detailed in Table 1. In this study, we also calculate total gravity offset index values and total gravity anomaly values 2 years before the earthquake in the southern Himalayan Orogenic Belt area. The results are shown in Figure 9. We find that the gravity offset index values and gravity anomaly values in the fault zone both appeared as phenomena of significant increases 6 months before the earthquake, and they reached a maximum during the period from 2013 to 2015. Figure 10 shows the spatial distribution of gravity offset index values in the southern Himalayan Orogenic Belt 6 months before the earthquake. As can be seen from Figure 10, significant gravity changes in the southern Himalayan Orogenic Belt occurred 3–6 months before the Nepal earthquake sequence, and the maximum K value was located approximately 50 km from the epicenter of the MW 7.8 Nepal earthquake, as shown in Figures 10b and 10c.
Figure 11 shows the spatiotemporal distribution of gravity anomalies whose K values are greater than 3. As shown in Figure 11b, a large area of anomalies appeared in the southern Himalayan Orogenic Belt region 6 months before the MW 7.8 Nepal earthquake. In particular, all these anomalous regions were located around six epicenters, and they occupied almost one third of the map area. The anomalies mentioned above indicate high stresses around the six epicenter locations. Thus, it can be inferred that the compression of the Indian Plate on the Tibetan Plateau was significantly enhanced from October 2014 to January 2015 and that the tectonic activities on both sides of the fault zone were extremely active. As can be seen in Figure 11f, multiple anomalies appeared along the southern Himalayan Orogenic Belt from March to May 2015, and the six epicenters of the Nepal earthquake sequence were all at locations with high anomalies. Considering that the gravity variation in GRACE data was positive during this period, it can be inferred that stress accumulation had been completed on the fault belt. However, the GRACE gravity variation of the fault zone and its southern region were negative when the earthquake happened, which suggests that the occurrence of the Nepal earthquake sequence was due to the fracture of the fault zone after the stress had reached a critical value. As shown in Figure 12, there were no gravity anomalies (the maximum K value was only 2.6) in the southern Himalayan Orogeny during the same period in nonearthquake years.
In sum, compared with a single earthquake, the seismic sequence can better reveal the process of tectonic activity in the fault zone before the earthquake. We find that before the earthquake, strong tectonic activities usually occur in the fault zone, thus bringing about the phenomenon of obvious stress concentration. Specifically, the distribution of the anomalies appeared along the fault zone, and the maximum anomaly occurred approximately 50 km from the epicenter half a year before the earthquake, suggesting it may be a kind of earthquake precursor signal.
Figure 13a displays all K value changes in the first half year of the Wenchuan earthquake, and Figures 13b and 13c show the spatiotemporal distribution of K values greater than 3 in this period. As can be seen in Figure 13a, the higher offset index K values on the Tibetan Plateau show a characteristic congregating distribution during the 6 months before the Wenchuan earthquake.
Figure 13b shows that from December 2007 to February 2008, almost half a year to 3 months before the Wenchuan earthquake, a wide range of NE-trending gravity anomalies first appeared in the southern Tibetan Plateau. In addition, a large area of abnormalities appeared in the northwest part of the plateau. One month later, the MW 7.2 Yutian earthquake occurred in this region. This was probably due to the materials on the plateau being blocked by the Tarim Basin in the process of migrating northward under the compression of the Indian Plate, causing a massive accumulation of material in the northwestern plateau. That accumulation allowed the stress to build up in this area, which eventually led to the occurrence of the Yutian earthquake.
At the same time, the theory of material eastward flow holds that since the collision between the Indian Plate and the Eurasian Plate, material on the Tibetan Plateau has flowed from west to east (Clark and Royden, 2000; Royden et al., 2008). Therefore, on the basis of this theory, it is not difficult to infer that when the material flows to the South China Block and is blocked by it, the direction of migration changes to the southeast, resulting in the accumulation of a large amount of material in the south of the South China Block. This inference is consistent with the large area of gravity anomalies that arose in the region during this time period. As the northwest boundary of the South China Block, the Longmenshan Fault Zone may have accumulated stress during the eastward movement of the crustal material of the Tibetan Plateau. This speculation can be confirmed by the occurrence of a large number of gravity anomalies in the Longmenshan Fault Zone 3 months later, as shown in Figure 13c.
Figure 13c shows the distribution of gravity anomalies in the study area during the 3 months before the Wenchuan earthquake. We find that the large range of gravity anomalies in the Himalayan Orogenic Belt is basically consistent with the trend of that mountain range, accounting for about one quarter of the total number of anomalies on the map. This result suggests that the plateau continues to experience compressional stress from the Indian Plate. The MW 7.2 Yutian earthquake occurred in the Altyn Fault Zone as extrusion stress was transferred to the inner block of the plateau. Xu XW et al. (2013) later confirmed that the crust rupture zone produced by the Yutian earthquake was a tensile shear rupture zone with both normal slip and left-lateral strike-slip components and that its seismogenic fault belongs to the NE-trending tensional tectonic zone at the southwest end of the Altyn Tagh Fault. Combined with the GPS data (Gan WJ et al., 2007), signs of the Bayan Har, Qaidam, and Qilian Blocks sliding eastward were clear during the period of this earthquake. Because the Longmenshan Fault Zone lies between these three blocks and the South China Block, a large amount of stress accumulated in the Longmenshan Fault Zone and its northwest side appeared blocked because of the strong obstruction by the South China Block, which provided the seismogenic conditions for the Wenchuan earthquake one month later.
Figure 14a shows all changes in K values in the first half year before the MW 7.8 Nepal earthquake, whereas Figures 14b and 14c display the spatiotemporal distribution of K values greater than 3 in this period. As shown in Figure 14a, the frequency of occurrence of high offset index K values between 6 and 3 months before the earthquake is significantly greater than that during the 3 months before the earthquake. Moreover, through further analysis of the high offset index K values, we find that between 6 and 3 months before the earthquake, several high anomalies greater than 4 were centered at the epicenter and distributed along the fault zone in the Himalayan Orogenic Belt. Simultaneously, a large area of gravity anomalies was also gathering on the south side of this belt.
The anomaly characterization method proposed in this study is based on the maximum shear strain. Although it can determine well where the anomaly exists, it cannot directly judge whether the abnormal region is due to stress extrusion or stress release. To this end, the variation in crustal mass density calculated based on GRACE data is analyzed here to help determine whether the abnormal area is a stress compressive zone or stress release zone. Generally, the crustal mass density increases in a compressive zone and decreases in a tension zone. Therefore, to explore the characteristics of mass density variation in the seismogenic region 2 years before the Nepal earthquake, we calculate the proportion, denoted as R, of the increased area of crustal mass density to the total area of the Indian Plate and the Tibetan Plateau, respectively.
As shown in Figure 15, the red and blue lines represent the change processes of R in the Tibetan Plateau and the Indian Plate, respectively, 2 years before the earthquake. During the 2 years before the earthquake, the crustal mass density of the Indian Plate was seen mainly to decrease, with the R value nearly always below 0.5, whereas between 6 and 3 months before the earthquake, the area of mass density increased continuously, until it covered 53% of the total area of the Indian Plate. Figure 16 is a visual presentation of the changes in crustal mass density of the study area in different time periods. As shown in Figure 16b, the mass density increased on the southwest side of the Himalayan Orogenic Belt from 6 to 3 months before the earthquake, indicating that this region was under a state of stress compression at this stage, which made it possible to trigger the subsequent Nepal earthquake. Gravity anomalies can also be observed near multiple faults on the Tibetan Plateau during this period. We speculated that the Nepal earthquake was triggered as a result of the Indian Plate pushing the Eurasian Plate where the Tibetan Plateau is located and that the force was transmitted to the interior of the plateau.
However, note that the area of gravity anomalies within the entire study area was significantly reduced during the 3 months before the earthquake, as shown in Figures 9 and 14c. This phenomenon can also be observed in Figure 15. As indicated by the blue line therein, the increasing area of crustal mass density within the Indian Plate decreased significantly during the 3 months before the earthquake. Instead, areas where mass density decreased were expanding during this period. Judging from the red curved segment in time period C in Figure 15, the material gradually flowed out from the Himalayan Orogenic Belt and south, then migrated into the interior of the Tibetan Plateau, under the continuous northward pushing effect of the Indian Plate to the Tibetan Plateau. As a result, the strain energy accumulated in the Himalayan Orogenic Belt and its southern regions was gradually released and the crustal mass density inside the Tibetan Plateau was increasing. According to the meta-instability theory proposed by Jin MA et al. (2012), the process of “stress accumulation and then release” indicates that the boundary between the Indian Plate and the Eurasian Plate had reached a state of instability during the 3 months before the earthquake. In the end, the Nepal earthquake broke out with a rapid release of stress on the fault zone. Moreover, Chen SY et al. (2020) utilized satellite remote sensing land surface temperature products to analyze the relationship between the stress change before the Nepal earthquake and the temperature response. They concluded that the observed coseismic temperature decrease was qualitatively in accordance with the coseismic stress drop on the upper side of the thrust fault. Considering that these changes occurred before the earthquake, they may even have served as precursors to the Nepal earthquake, further confirming the speculation above.
In this study, no filtering operation is carried out in the process of calculating the maximum shear strain of gravity, which eventually avoids the loss of higher order information in GRACE data. In addition, the maximum shear strain of gravity is calculated based on the gravity gradient, which has a good suppression effect on the striping noise (Wang L et al., 2012; Li J et al., 2015; Dai CL et al., 2016). In other words, the maximum shear strain of gravity can effectively suppress noise without losing gravity information. To further verify the effectiveness of the proposed method, four commonly used GRACE data processing methods are used as comparisons to comprehensively analyze the preseismic gravity anomaly extraction results obtained from multiple earthquake cases. These comparison methods include three Gaussian filtering methods with filter radii of 250, 300, and 500 km, respectively, and one postprocessing method combining Gaussian filtering (with a filter radius of 300 km) and a correlated-error filter.
The statistical results in Table 2 show that the proposed method can extract pre-earthquake gravity anomaly information for all the earthquake cases with a magnitude of 7.5 and above, whereas the highest successful extraction rate among the comparison methods is only 50%. For earthquake cases with a magnitude of 6 and above, the highest successful extraction rate among the existing methods is only 38.5%, whereas the proposed method is as high as 61.5%, which is 1.6 times that of the existing methods. This result fully demonstrates that the method of extracting gravity anomalies before an earthquake without filtering can greatly improve the sensitivity of gravity anomaly information. In addition, it can be concluded that Gaussian filtering with a filtering radius of 500 km will cause the loss of seismic signals in GRACE data, which is basically consistent with the results of Sun WK et al. (2011) and Zhao Q et al. (2015). In conclusion, the method presented in this article can extract the pre-earthquake gravity anomaly information well for large earthquakes with a magnitude of 7.5 or above.
ID | Event | K values | Preseismic gravity values based on existing methods | ||||||||
Gaussian filter (R = 250 km) |
Gaussian filter (R = 300 km) |
Gaussian filter (R = 300 km) + correlated-error filter |
Gaussian filter (R = 500 km) |
||||||||
A1 | T2 | A | T | A | T | A | T | A | T | ||
1 | Wenchuan MW 7.9 | Y | 1−3 | N | N | N | N | ||||
2 | Sichuan MW 6.1 | Y | From −1 to 1 | N | N | N | N | ||||
3 | Sichuan–Gansu MW 6.1 | Y | 1−3 | N | N | N | N | ||||
4 | Qingchuan MW 6.0 | Y | 5−7 | N | N | N | N | ||||
5 | Xizang MW 6.1 | N | N | N | N | N | |||||
6 | Yutian MW 7.2 | N | N | N | N | N | |||||
7 | Yutian MW 6.2 | N | N | N | N | N | |||||
8 | Yutian MW 6.9 | Y | 3−5 | N | N | N | N | ||||
9 | Nepal MW 7.8 | Y | 3−6 | Y | 3 | Y | 3 | Y | 3 | N | |
10 | Nepal MW 6.6 | Y | 3−6 | Y | 3 | Y | 3 | Y | 3 | N | |
11 | Nepal MW 6.1 | Y | 9−11 | Y | 3 | Y | 3 | Y | 3 | N | |
12 | Nepal MW 7.3 | N | N | Y | 3 | Y | 3 | N | |||
13 | Nepal MW 6.3 | N | N | Y | 3 | Y | 3 | N | |||
8 Y 5 N 100% Y (for MW ≥ 7.5) 61.5% Y (for MW ≥ 6) |
3 Y 10 N 50% Y (for MW ≥ 7.5) 23.1% Y (for MW ≥ 6) |
5 Y 8 N 50% Y (for MW ≥ 7.5) 38.5% Y (for MW ≥ 6) |
5 Y 8 N 50% Y (for MW ≥ 7.5) 38.5% Y (for MW ≥ 6) |
0 Y 13 N 0% Y (for MW ≥ 7.5) 0% Y (for MW ≥ 6) |
|||||||
1A, anomaly. The presence (Y) or absence (N) of gravity anomalies in the time series is determined. 2T, time. The unit of time is a month, and the occurrence time of the earthquake is the origin time. A positive number indicates before the earthquake, and a negative number indicates after the earthquake. For example, from −1 to 1 indicates one month before and one month after the earthquake. |
Some scholars believe that tectonic plate movements and the interactions between plates are the main factors contributing to seismic development and occurrence (Keilis, 1990; Zhang PZ et al., 2003; Zhang GM, 2005; Ismail-Zadeh et al., 2007). And according to Equation (11), when the magnitude is 8, the radius of the seismogenic zone is more than 2700 km, whereas when the magnitude is 9, the radius of the seismogenic zone can increase to more than 7000 km. Therefore, the gravity variation in the Tibetan Plateau and the Indian Plate, 3000 km or less away from the Sunda Subplate, can be inferred to have been abnormal during the 6 months before the famous MW 9.1 Sumatra earthquake occurred on the Sunda Subplate in 2004. To confirm this hypothesis, we take the occurrence of large earthquakes as the observed time node and investigate the temporal synchronization of pre-earthquake gravity anomalies between the Tibetan Plateau and its adjacent plates, such as the Indian Plate and the Sunda Subplate, as shown in Figure 17. For this reason, this study draws the variation curve of total K values for all pixels in the regions of the Tibetan Plateau, Indian Plate, and Sunda Subplate, respectively, from 2003 to 2016. The results are illustrated in Figure 18. As marked by the red circles, 6 months before the Nepal earthquake occurred at the boundary between the Tibetan Plateau and the Indian Plate, K values on the Tibetan Plateau, Indian Plate, and Sunda Subplate all underwent a rapid increase simultaneously and soon reached their maximum values. Furthermore, 3 months prior to the MW 9.1 Sumatra earthquake that occurred in the Sunda Subplate, synchronous and remarkable increases in K values were observed in the Tibetan Plateau, Indian Plate, and Sunda Subplate, and all reached their maxima during the period from 2003 to 2014, as expressed by the blue circle in Figure 18.
In sum, it can be supposed that a strong correlation of tectonic activities exists between adjacent plates, based on the fact that gravity anomalies on adjacent plates were synchronized before large earthquakes. Moreover, it confirms that the Nepal and Sumatra earthquakes may be related to the enhanced interactions between the Tibetan Plateau, Indian Plate, and Sunda Subplate.
The GRACE system is widely used in earthquake monitoring because of its great ability to detect changes in the seismic gravitational field. In this study, we propose a method of extracting the pre-earthquake gravity anomaly information based on the maximum shear strain. Exploratory experiments are carried out on the Tibetan Plateau and the surrounding areas, located among the Pacific, Indian, and Eurasian Plates, with the highest altitude, most complex topography, and most frequent strong earthquakes. Ultimately, the spatiotemporal characteristics of gravity anomalies before the Wenchuan earthquake and Nepal earthquake are thoroughly analyzed, and the main conclusions can be summarized as follows:
(1) At the observation scale of the tectonic fault zone, we find that gravity anomalies of large areas occurred nearby the epicenter many times about half a year before the earthquake and that these anomalies were distributed along the fault zone. At the observation scale of the plate, we find that when an earthquake occurred on the Tibetan Plateau, a large number of gravity anomalies also appeared at the boundary of the Tibetan Plateau and the Indian Plate. Conversely, unlike the gravity anomalies of large areas appearing in the preseismic period, no obvious gravity anomalies appeared during the nonseismic period, no matter whether at the plate or fault zone scale. This finding indicates that the tectonic activities in intraplate and interplate regions before an earthquake were more active than those in the nonseismic period.
(2) Because the maximum shear strain is calculated by using the second-order gradient of the disturbance potential, the proposed method can effectively suppress the stripe noise from GRACE data, thus greatly increasing the sensitivity of detecting pre-seismic gravity anomalies. Therefore, it becomes possible to detect the pre-seismic gravity anomalies of earthquakes with magnitudes of less than 8, such as the Wenchuan earthquake and the Nepal earthquake. In general, this study provides new ideas for applying gravity satellite data to earthquake monitoring research.
We gratefully appreciate the editor and anonymous reviewers for their efforts and constructive comments, which have greatly improved the technical quality and presentation of this study.
Allen, R. M., and Melgar, D. (2019). Earthquake early warning: Advances, scientific challenges, and societal needs. Annu. Rev. Earth Planet. Sci., 47, 361–388. https://doi.org/10.1146/annurev-earth-053018-060457
|
Barnes, D. F. (1966). Gravity changes during the Alaska earthquake. J. Geophys. Res., 71(2), 451–456. https://doi.org/10.1029/jz071i002p00451
|
Brozzetti, F., Mondini, A. C., Pauselli, C., Mancinelli, P., Cirillo, D., Guzzetti, F., and Lavecchia, G. (2020). Mainshock anticipated by intra-sequence ground deformations: Insights from multiscale field and SAR interferometric measurements. Geosciences, 10(5), 186. https://doi.org/10.3390/geosciences10050186
|
Burg, J. P., Davy, P., Nievergelt, P., Oberli, F., Seward, D., Diao, Z. Z., and Meier, M. (1997). Exhumation during crustal folding in the Namche-Barwa syntaxis. Terra Nova, 9(2), 53–56. https://doi.org/10.1111/j.1365-3121.1997.tb00001.x
|
Chao, B. F., and Liau, J. R. (2019). Gravity changes due to large earthquakes detected in GRACE satellite data via empirical orthogonal function analysis. J. Geophys. Res.: Solid Earth, 124(3), 3024–3035. https://doi.org/10.1029/2018JB016862
|
Chen, S., Liu, M., Xing, L. L., Xu, W. M., Wang, W. X., Zhu, Y. Q., and Li, H. (2016). Gravity increase before the 2015 MW 7.8 Nepal earthquake. Geophys. Res. Lett., 43(1), 111–117. https://doi.org/10.1002/2015GL066595
|
Chen, S. Y., Liu, P. X., Feng, T., Wang, D., Jiao, Z. H., Chen, L. C., Xu, Z. X., and Zhang, G. Z. (2020). Exploring changes in land surface temperature possibly associated with earthquake: Case of the April 2015 Nepal MW 7.9 earthquake. Entropy, 22(4), 377. https://doi.org/10.3390/e22040377
|
Cheng, M. K., and Tapley, B. D. (2004). Variations in the Earth’s oblateness during the past 28 years. J. Geophys. Res.: Solid Earth, 109(B9), B09402. https://doi.org/10.1029/2004JB003028
|
Clark, M. K., and Royden, L. H. (2000). Topographic ooze: Building the eastern margin of Tibet by lower crustal flow. Geology, 28(8), 703. https://doi.org/10.1130/0091-7613(2000)28<703:TOBTEM>2.0.CO;2 doi: 10.1130/0091-7613(2000)28<703:TOBTEM>2.0.CO;2
|
Cui, J., Chen, W., Li, P., Zhang, X., and Li, L. (2003). Thrust propagation in the Aqqikkol Lake Area, the East Kunlun Mountains, Northwestern China. Acta Geologica Sinica-English Edition, 77(4), 468–478. https://doi.org/10.1111/j.1755-6724.2003.tb00127.x
|
D’Alessandro, A., D’Anna, R., Greco, L., Passafiume, G., Scudero, S., Speciale, S., and Vitale, G. (2018). Monitoring earthquake through MEMS sensors (MEMS project) in the town of Acireale (Italy). In 5th IEEE International Symposium on Inertial Sensors and Systems (pp. 1–4). Lake Como, Italy: IEEE. https://doi.org/10.1109/ISISS.2018.8358143
|
Dai, C. L., Shum, C. K., Guo, J. Y., Shang, K., Tapley, B., and Wang, R. J. (2016). Improved source parameter constraints for five undersea earthquakes from north component of GRACE gravity and gravity gradient change measurements. Earth Planet. Sci. Lett., 443, 118–128. https://doi.org/10.1016/j.jpgl.2016.03.025
|
Dermanis, A., and Livieratos, E. (1983). Applications of deformation analysis in geodesy and geodynamics. Rev. Geophys., 21(1), 41–50. https://doi.org/10.1029/RG021i001p00041
|
Dirks, P. H. G. M., Wilson, C. J. L., Chen, S., Luo, Z. L., and Liu, S. (1994). Tectonic evolution of the NE margin of the Tibetan Plateau; evidence from the central Longmen Mountains, Sichuan Province, China. J. Southeast Asian Earth Sci., 9(1-2), 181–192. https://doi.org/10.1016/0743-9547(94)90074-4
|
Dobrovolsky, I. P., Zubkov, S. I., and Miachkin, V. I. (1979). Estimation of the size of earthquake preparation zones. Pure Appl. Geophys., 117(5), 1025–1044. https://doi.org/10.1007/BF00876083
|
Duan, H. R., Guo, J. G., Chen, L. K., Jiao, J. S., and Jian, H. T. (2022). Vertical crustal deformation velocity and its influencing factors over the Qinghai–Tibet Plateau based on satellite gravity data. Earth Planet. Phys., 6(4), 366–377 https://doi.org/10.26464/epp2022034
|
Eshagh, M. (2008). Non-singular expressions for the vector and the gradient tensor of gravitation in a geocentric spherical frame. Comput. Geosci., 34(12), 1762–1768. https://doi.org/10.1016/j.cageo.2008.02.022
|
Eshagh, M., and Abdollahzadeh, M. (2012). Software for generating gravity gradients using a geopotential model based on an irregular semivectorization algorithm. Comput. Geosci., 39, 152–160. https://doi.org/10.1016/j.cageo.2011.06.003
|
Fan, X. M., Scaringi, G., Korup, O., West, A. J., Van Westen, C. J., Tanyas, H., Hovius, N., Hales, T. C., Jibson, R. W., … Huang, R. Q. (2019). Earthquake-induced chains of geologic hazards: Patterns, mechanisms, and impacts. Rev. Geophys., 57(2), 421–503. https://doi.org/10.1029/2018RG000626
|
Fatolazadeh, F., Goïta, K., and Azar, R. J. (2020). Determination of earthquake epicentres based upon invariant quantities of GRACE strain gravity tensors. Sci. Rep., 10, 7636. https://doi.org/10.1038/s41598-020-64560-w
|
Gan, W. J., Zhang, P. Z., Shen, Z. K., Niu, Z. J., Wang, M., Wan, Y. G., Zhou, D. M., and Cheng, J. (2007). Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements. J. Geophys. Res.: Solid Earth, 112(B8), B08416. https://doi.org/10.1029/2005JB004120
|
Gido, N. A. A., Amin, H., Bagherbandi, M., and Nilfouroushan, F. (2020). Satellite monitoring of mass changes and ground subsidence in Sudan’s oil fields using GRACE and Sentinel-1 data. Remote Sens., 12(11), 1792. https://doi.org/10.3390/rs12111792
|
Han, S. C., Shum, C. K., Bevis, M., Ji, C., and Kuo, C. Y. (2006). Crustal dilatation observed by GRACE after the 2004 Sumatra-Andaman earthquake. Science, 313(5787), 658–662. https://doi.org/10.1126/science.1128661
|
Ince, E. S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F., and Schuh, H. (2019). ICGEM–15 years of successful collection and distribution of global gravitational models, associated services, and future plans. Earth Syst. Sci. Data, 11(2), 647–674. https://doi.org/10.5194/essd-11-647-2019
|
Ismail-Zadeh, A., Le Mouël, J. L., Soloviev, A., Tapponnier, P., and Vorovieva, I. (2007). Numerical modeling of crustal block-and-fault dynamics, earthquakes and slip rates in the Tibet-Himalayan region. Earth Planet. Sci. Lett., 258(3-4), 465–485. https://doi.org/10.1016/j.jpgl.2007.04.006
|
Jin, M. A., Sherman, S. I., and Guo, Y. S. (2012). Identification of meta-instable stress state based on experimental study of evolution of the temperature field during stick-slip instability on a 5° bending fault. Sci. China Earth Sci., 55(6), 869–881. https://doi.org/10.1007/s11430-012-4423-2
|
Keilis-Borok, V. I. (1990). The lithosphere of the Earth as a nonlinear system with implications for earthquake prediction. Rev. Geophys., 28(1), 19–34. https://doi.org/10.1029/RG028i001p00019
|
Li, J., and Shen, W. B. (2015). Monthly GRACE detection of coseismic gravity change associated with 2011 Tohoku-Oki earthquake using northern gradient approach. Earth, Planets Sp., 67(1), 29. https://doi.org/10.1186/s40623-015-0188-0
|
Li, L. M., Wen, Z. Z., and Wang Z. S. (2016). Outlier detection and correction during the process of groundwater lever monitoring base on Pauta criterion with self-learning and smooth processing. In 16th Asian Simulation Conference, SCS Autumn Simulation Multi-Conference (pp. 497–503). Beijing, China: Springer. https://doi.org/10.1007/978-981-10-2663-8_51
|
Liu, C., and Sun, W. K. (2023). GRACE time-variable gravity and it sapplication to geoscience: Quantitative analysis of relevant literature. Earth Planet. Phys., 7(2), 295–309. https://doi.org/10.26464/epp2023022 https://doi.org/10.26464/epp2023022
|
Mahmood, I. (2021). Anomalous variations of air temperature prior to earthquakes. Geocarto Int., 36(12), 1396–1408. https://doi.org/10.1080/10106049.2019.1648565
|
Marchetti, D., De Santis, A., Shen, X. H., Campuzano, S. A., Perrone, L., Piscini, A., Di Giovambattista, R., Jin, S. G., Ippolito, A., … Huang, J. P. (2020). Possible lithosphere-atmosphere-ionosphere coupling effects prior to the 2018 MW = 7.5 Indonesia earthquake from seismic, atmospheric and ionospheric data. J. Asian Earth Sci., 188, 104097. https://doi.org/10.1016/j.jseaes.2019.104097
|
Martinelli, G. (2020). Previous, current, and future trends in research into earthquake precursors in geofluids. Geosciences, 10(5), 189. https://doi.org/10.3390/geosciences10050189
|
Massonnet, D., Feigl, K., Rossi, M., and Adragna, F. (1994). Radar interferometric mapping of deformation in the year after the Landers earthquake. Nature, 369(6477), 227–230. https://doi.org/10.1038/369227a0
|
Matsuo, K., and Heki, K. (2011). Coseismic gravity changes of the 2011 Tohoku-Oki earthquake from satellite gravimetry. Geophys. Res. Lett., 38(7), L00G12. https://doi.org/10.1029/2011GL049018
|
Mikhailov, V. O., Timoshkina, E. P., Smirnov, V. B., Khairetdinov, S. A., and Dmitriev, P. N. (2020). On the origin of postseismic deformation processes in the region of the Maule, Chile earthquake of February 27, 2010. Izv. Phys. Solid Earth, 56(6), 762–771. https://doi.org/10.1134/S106935132006004X
|
Molnar, P., and Lyon-Caent, H. (1989). Fault plane solutions of earthquakes and active tectonics of the Tibetan Plateau and its margins. Geophys. J. Int., 99(1), 123–154. https://doi.org/10.1111/j.1365-246X.1989.tb02020.x
|
Panet, I., Pollitz, F., Mikhailov, V., Diament, M., Banerjee, P., and Grijalva, K. (2010). Upper mantle rheology from GRACE and GPS postseismic deformation after the 2004 Sumatra-Andaman earthquake. Geochem. Geophys. Geosyst., 11(6), Q06008. https://doi.org/10.1029/2009GC002905
|
Panet, I., Bonvalot, S., Narteau, C., Remy, D., and Lemoine, J. M. (2018). Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data. Nat. Geosci., 11(5), 367–373. https://doi.org/10.1038/s41561-018-0099-3
|
Peidou, A., and Pagiatakis, S. (2019). Gravity gradiometry with GRACE space missions: New opportunities for the geosciences. J. Geophys. Res.: Solid Earth, 124(8), 9130–9147. https://doi.org/10.1029/2018JB016382
|
Peidou, A., and Pagiatakis, S. (2020). Stripe mystery in GRACE geopotential models revealed. Geophys. Res. Lett., 47(4), e2019GL085497. https://doi.org/10.1029/2019GL085497
|
Petrovskaya, M. S., and Vershkov, A. N. (2008). Development of the second-order derivatives of the Earth’s potential in the local north-oriented reference frame in orthogonal series of modified spherical harmonics. J. Geod., 82(12), 929–944. https://doi.org/10.1007/s00190-008-0223-z
|
Qu, W., Lu, Z., Zhang, Q., Hao, M., Wang, Q. L., Qu, F. F., and Zhu, W. (2018). Present-day crustal deformation characteristics of the southeastern Tibetan Plateau and surrounding areas by using GPS analysis. J. Asian Earth Sci., 163, 22–31. https://doi.org/10.1016/j.jseaes.2018.05.021
|
Rao, W. L., and Sun, W. K. (2022). Runoff variations in the Yangtze River Basin and sub-basins based on GRACE, hydrological models, and in-situ data. Earth Planet. Phys., 6(3), 228–240. https://doi.org/10.26464/epp2022021 https://doi.org/10.26464/epp2022021
|
Rodell, M., Houser, P. R., Jambor, U., Gottschalck, J., Mitchell, K., Meng, C. J., Arsenault, K., Cosgrove, B., Radakovich, J., .. Toll, D. (2004). The global land data assimilation system. Bulletin of the American Meteorological Society, 85(3), 381–394. https://doi.org/10.1175/BAMS-85-3-381
|
Royden, L. H., Burchfiel, B. C., and Van Der Hilst, R. D. (2008). The geological evolution of the Tibetan Plateau. Science, 321, 1054–1058. https://doi.org/10.1126/science.1155371
|
Saç, M. M., Harmanşah, C., Camgöz, B., and Sözbilir, H. (2011). Radon monitoring as the earthquake precursor in fault line in Western Turkey. Ekoloji, 20(79), 93–98. https://doi.org/10.5053/ekoloji.2011.7912
|
Sansò, F., and Barzaghi, R. (2020). Foreword: Earth’s gravity field and Earth sciences. Rend. Lincei. Sci. Fis. Natur., 31(1), 1–2. https://doi.org/10.1007/s12210-020-00961-3
|
Sasmal, S., Chowdhury, S., Kundu, S., Politis, D. Z., Potirakis, S. M., Balasis, G., Hayakawa, M., and Chakrabarti, S. K. (2021). Pre-seismic irregularities during the 2020 Samos (Greece) earthquake (M = 6.9) as investigated from multi-parameter approach by ground and space-based techniques. Atmosphere (Basel), 12(8), 1059. https://doi.org/10.3390/atmos12081059
|
Seheult, A. H., Green, P. J., Rousseeuw, P. J., and Leroy, A. M. (1989). Robust regression and outlier detection. J. Roy. Stat. Soc. Ser. A Stat. Soc., 152(1), 133. https://doi.org/10.2307/2982847
|
Song, D. M., Xie, R. H., Zang, L., Yin, J. Y., Qin, K., Shan, X. J., Cui, J. Y., and Wang, B. (2018). A new algorithm for the characterization of thermal infrared anomalies in tectonic activities. Remote Sens., 10(12), 1941. https://doi.org/10.3390/rs10121941
|
Sun, H. P., Xu, J. Q., and Cui, X. M. (2017). Research progress of the gravity field application in Earth’s geodynamics and interior structure. Acta Geod. Cartogr. Sin. (in Chinese), 46(10), 1290–1299. https://doi.org/10.11947/J.AGCS.2017.20170298
|
Sun, W., and Okubo, S. (2004). Coseismic deformations detectable by satellite gravity missions: A case study of Alaska (1964, 2002) and Hokkaido (2003) earthquakes in the spectral domain. J. Geophys. Res.: Solid Earth, 109(B4), B04405. https://doi.org/10.1029/2003jb002554
|
Sun, W. K., Hasegawa, T., Zhang, X. L., Fukuda, Y., Shum, C. K., and Wang, L. (2011). Effects of Gaussian filter in processing GRACE data: Gravity rate of change at Lhasa, southern Tibet. Sci. China Earth Sci., 54(9), 1378–1385. https://doi.org/10.1007/s11430-011-4233-y
|
Swenson, S., Chambers, D., and Wahr, J. (2008). Estimating geocenter variations from a combination of GRACE and ocean model output. J. Geophys. Res.: Solid Earth, 113(B8), B08410. https://doi.org/10.1029/2007JB005338
|
Tan, C. X., Zhang, P., Lu, S. L., Zhu, J. Z., Feng, C. J., Qin, X. H., and Meng, J. (2019). Significance and role of in-situ crustal stress measuring and real-time monitoring in earthquake prediction research. J. Geomech. (in Chinese), 25(5), 866–876. https://doi.org/10.12090/j.issn.1006-6616.2019.25.05.071
|
Tang, H., and Sun, W. K. (2018). Closed-form expressions of seismic deformation in a homogeneous Maxwell earth model. J. Geophys. Res.: Solid Earth, 123(7), 6033–6051. https://doi.org/10.1029/2018JB015594
|
Tang, H., and Sun, W. (2021). Analytical equations for an infinite series involving low-order associated Legendre functions in geoscience. J. Geod., 95(7), 86. https://doi.org/10.1007/s00190-021-01527-3
|
Tang, L., Li, J., Chen, J. L., Wang, S. Y., Wang, R., and Hu, X. G. (2020). Seismic impact of large earthquakes on estimating global mean ocean mass change from GRACE. Remote Sens., 12(6), 935. https://doi.org/10.3390/rs12060935
|
Taylor, M., and Yin, A. (2009). Active structures of the Himalayan-Tibetan orogen and their relationships to earthquake distribution, contemporary strain field, and Cenozoic volcanism. Geosphere, 5(3), 199–214. https://doi.org/10.1130/GES00217.1
|
Tramutoli, V. (1998). Robust AVHRR techniques (RAT) for environmental monitoring: Theory and applications. In Proceedings of SPIE 3496, Earth Surface Remote Sensing II. Barcelona: SPIE. https://doi.org/10.1117/12.332714
|
Velicogna, I., Mohajerani, Y., Geruo, A. G., Landerer, F., Mouginot, J., Noel, B., Rignot, E., Sutterley, T., Van Den broeke, M., … Wiese, D. (2020). Continuity of ice sheet mass loss in Greenland and Antarctica from the GRACE and GRACE follow-on missions. Geophys. Res. Lett., 47(8), e2020GL087291. https://doi.org/10.1029/2020GL087291
|
Wahr, J., Molenaar, M., and Bryan, F. (1998). Time variability of the Earth’s gravity field: Hydrological and oceanic effects and their possible detection using GRACE. J. Geophys. Res.: Solid Earth, 103(B12), 30205–30229. https://doi.org/10.1029/98jb02844
|
Wang, L., Shum, C. K., and Jekeli, C. (2012). Gravitational gradient changes following the 2004 December 26 Sumatra-Andaman Earthquake inferred from GRACE. Geophys. J. Int., 191(3), 1109–1118. https://doi.org/10.1111/j.1365-246X.2012.05674.x
|
Wang, L. H., Chen, S., Zhuang, J. C., Zhang, B., Shi, W., Yang, J. L., and Xu, W. M. (2023). Gravity field changes reveal deep mass transfer before and after the 2013 Lushan earthquake. Commun. Earth Environ., 4(1), 194. https://doi.org/10.1038/s43247-023-00860-z
|
Wang, X. L., Luo, Z. C., Zhong, B., Wu, Y. H., Huang, Z. K., Zhou, H., and Li, Q. (2019). Separation and recovery of geophysical signals based on the Kalman filter with GRACE gravity data. Remote Sens., 11(4), 393. https://doi.org/10.3390/rs11040393
|
Xu, X. W., Tan, X. B., Yu, G. H., Wu, G. D., Fang, W., Chen, J. B., Song, H. P., and Shen, J. (2013). Normal- and oblique-slip of the 2008 Yutian earthquake: Evidence for eastward block motion, northern Tibetan Plateau. Tectonophysics, 584, 152–165. https://doi.org/10.1016/j.tecto.2012.08.007
|
Zhang, G. M., Ma, H. S., Wang, H., and Wang, X. L. (2005). Boundaries between active-tectonic blocks and strong earthquakes in China mainland. Chin. J. Geophys., 48(3), 662–671. https://doi.org/10.1002/cjg2.699
|
Zhang, P. Z., Deng, Q. D., Zhang, G. M., Ma, J., Gan, W. J., Min, W., Mao, F. Y., and Wang, Q. (2003). Active tectonic blocks and strong earthquakes in the continent of China. Sci. China, Ser. D Earth Sci., 46(S2), 13–24. https://doi.org/10.1360/03dz0002
|
Zhang, Y., Wu, Y. L., Yan, J. G., Wang, H. R., Rodriguez, J. A. P., and Qiu, Y. (2018). 3D inversion of full gravity gradient tensor data in spherical coordinate system using local north-oriented frame. Earth, Planets Sp., 70(1), 58. https://doi.org/10.1186/s40623-018-0825-5
|
Zhao, Q., Wu, Y. L., and Wu, W. W. (2015). Effectiveness of empirical orthogonal function used in decorrelation of GRACE time-variable gravity field. Geod. Geodynam., 6(5), 324–332. https://doi.org/10.1016/j.geog.2015.05.003
|
Zhao, X. W., Pan, S., Sun, Z. C., Guo, H. D., Zhang, L., and Feng, K. (2021). Advances of satellite remote sensing technology in earthquake prediction. Nat. Hazards Rev., 22(1), 03120001. https://doi.org/10.1061/(asce)nh.1527-6996.0000419
|
Zhu, Y. Q., Yang, X., Liu, F., Zhao, Y. F., Wei, S. C., and Zhang, G. Q. (2023). Progress and prospect of the time-varying gravity in earthquake prediction in the Chinese Mainland. Front. Earth Sci., 11, 1124573. https://doi.org/10.3389/feart.2023.1124573
|
Zieger, T., and Ritter, J. R. R. (2018). Influence of wind turbines on seismic stations in the upper Rhine graben, SW Germany. J. Seismol., 22(1), 105–122. https://doi.org/10.1007/s10950-017-9694-9
|
Zotin, A., Simonov, K., Matsulev, A., and Kurako, M. (2019). Evaluation of gravitational anomalies in the areas of strongest earthquakes based on GRACE satellite measurements. Proc. Comput. Sci., 159, 1642–1651. https://doi.org/10.1016/j.procs.2019.09.334
|
Seismogenic fault | Data | Epicenter location | MW | Depth (km) | |
Longitude (°E) | Latitude (°N) | ||||
Longmenshan | 2008 May 12 | 103.322 | 31.002 | 7.9 | 19.0 |
2008 May 12 | 103.618 | 31.214 | 6.1 | 10.0 | |
2008 May 25 | 105.423 | 32.56 | 6.1 | 18.0 | |
2008 August 12 | 105.494 | 32.756 | 6.0 | 6.0 | |
Himalayan Orogenic Belt | 2015 April 25 | 84.731 | 28.231 | 7.8 | 8.2 |
2015 April 25 | 85.540 | 27.629 | 6.1 | 10.0 | |
2015 April 25 | 84.822 | 28.224 | 6.6 | 10.0 | |
2015 April 26 | 86.017 | 27.771 | 6.7 | 22.9 | |
2015 May 12 | 86.066 | 27.809 | 7.3 | 15.0 | |
2015 May 12 | 86.162 | 27.625 | 6.3 | 15.0 |
ID | Event | K values | Preseismic gravity values based on existing methods | ||||||||
Gaussian filter (R = 250 km) |
Gaussian filter (R = 300 km) |
Gaussian filter (R = 300 km) + correlated-error filter |
Gaussian filter (R = 500 km) |
||||||||
A1 | T2 | A | T | A | T | A | T | A | T | ||
1 | Wenchuan MW 7.9 | Y | 1−3 | N | N | N | N | ||||
2 | Sichuan MW 6.1 | Y | From −1 to 1 | N | N | N | N | ||||
3 | Sichuan–Gansu MW 6.1 | Y | 1−3 | N | N | N | N | ||||
4 | Qingchuan MW 6.0 | Y | 5−7 | N | N | N | N | ||||
5 | Xizang MW 6.1 | N | N | N | N | N | |||||
6 | Yutian MW 7.2 | N | N | N | N | N | |||||
7 | Yutian MW 6.2 | N | N | N | N | N | |||||
8 | Yutian MW 6.9 | Y | 3−5 | N | N | N | N | ||||
9 | Nepal MW 7.8 | Y | 3−6 | Y | 3 | Y | 3 | Y | 3 | N | |
10 | Nepal MW 6.6 | Y | 3−6 | Y | 3 | Y | 3 | Y | 3 | N | |
11 | Nepal MW 6.1 | Y | 9−11 | Y | 3 | Y | 3 | Y | 3 | N | |
12 | Nepal MW 7.3 | N | N | Y | 3 | Y | 3 | N | |||
13 | Nepal MW 6.3 | N | N | Y | 3 | Y | 3 | N | |||
8 Y 5 N 100% Y (for MW ≥ 7.5) 61.5% Y (for MW ≥ 6) |
3 Y 10 N 50% Y (for MW ≥ 7.5) 23.1% Y (for MW ≥ 6) |
5 Y 8 N 50% Y (for MW ≥ 7.5) 38.5% Y (for MW ≥ 6) |
5 Y 8 N 50% Y (for MW ≥ 7.5) 38.5% Y (for MW ≥ 6) |
0 Y 13 N 0% Y (for MW ≥ 7.5) 0% Y (for MW ≥ 6) |
|||||||
1A, anomaly. The presence (Y) or absence (N) of gravity anomalies in the time series is determined. 2T, time. The unit of time is a month, and the occurrence time of the earthquake is the origin time. A positive number indicates before the earthquake, and a negative number indicates after the earthquake. For example, from −1 to 1 indicates one month before and one month after the earthquake. |