Far‐field coseismic gravity changes related to the 2015 MW7.8 Nepal (Gorkha) earthquake observed by superconducting gravimeters in mainland China

Using data from five SGs at four stations in mainland China, obvious permanent gravity changes caused by the 2015 MW 7.8 Nepal (Gorkha) earthquake were detected. We analyzed the gravity effects from ground vertical deformation (VD) using co‐site continuous GPS (cGPS) data collocated at the Lijiang and the Wuhan station, and hydrological effects using GLDAS models and groundwater level records. After removing these effects, SG observations before and after the earthquake revealed obvious permanent gravity changes: −3.0 μGal, 7.3 μGal and 8.0 μGal at Lhasa, Lijiang and Wuhan station, respectively. We found that the gravity changes cannot be explained by the results of dislocation theory.


Introduction
Ground-based continuous gravity measurements reflect a combination of the transport and exchange of material and surface and interior deformation of the Earth, which are influenced by various environmental perturbations and geodynamic processes. Gravity changes recorded by a gravimeter consist of two factors: the contribution of internal mass redistribution and the effect of vertical displacement of Earth's surface. Earthquakes in general are instantaneous events, which redistribute local mass. The larger the earthquake, the greater the redistribution and relocation of mass around the earthquake epicenter. During an earthquake, abrupt frictional sliding along a part of the fault zone causes coseismic deformations and mass distributions (Van Camp et al., 2017). The amplitude of the observed coseismic gravity change depends on the magnitude of the earthquake and the distance between the epicenter and the site of gravity observation (Prasad et al., 2017). The high sensitivity, stability, and broad frequency range of dynamic linear responses of SG allow us to observe changes in gravity acceleration on the order of 1 μGal (1× 10 −8 ms −2 ) or even smaller (Imanishi et al., 2009). Tanaka et al. (2001) observed earthquake-induced gravity changes with an absolute gravimeter for the first time, just one day before and only 7 days after a M6.1 earthquake. The absolute gravity change was −6 μGal, significantly larger than the observational error of ~1 μGal. Imanishi et al. (2004) observed coseismic gravity changes by using an array of SGs located at epicentral distances of 3.4°, 6.9° and 9.4° from the 2003 M W 8.0 Tokachi-Oki earthquake. They observed gravity changes at the sub-microgal level, and demonstrated the ability of gravity measurements to infer the nature of dislocations in the earthquake source region. Both results are in good agreement with the calculated results of dislocation models. Kim et al. (2009) pointed out that the analysis period should not be too large, in order to minimize interaction with gravity variations from other sources. Nawa et al. (2009) detected gravity changes caused by the 2004 earthquakes off the Kii peninsula (M W 7.1, M W 7.4) using SG data from the Inuyama station, Japan. Prasad et al. (2017) observed a few microgal coseismic temporal gravity changes for the 2011 Koyna M W 4.8 earthquake in gPhone data, however, coseismic gravity changes calculated from dislocation theory (~0.1 μGal) were not in agreement. This implies that all the complex fault structures and perhaps weak fracture zones in the region below Koyna and Warna reservoirs were responding quickly to the accumulated stress, with reservoir loading and mass redistribution around the region resulting in a huge change in the gravity anomaly. Gravity changes in the above studies were smaller than 1 μGal with epicentral distances less than 100 km, and the earthquakes studied were typical thrust events occurring along oceanic trenches.
Estimation of the regional to global-scale gravity changes induced by a great earthquake requires consideration of the curvature and radial inhomogeneity of the Earth because these unusually large fault areas cover several thousands of kilometers. Sun WK andOkubo (1993, 1998) succeeded in obtaining the theoretical coseismic gravity changes for a spherically symmetric, non-rotating, perfectly elastic, and isotropic (SNREI) Earth with self-gravitation. Further, the GRACE (Gravity Recovery and Climate Experiment) satellite mission detected near-field significant gravity changes during the 2004 M W 9.1 Indian Ocean earthquake off the coast of Sumatra, Indonesia (Han et al., 2006), the 2010 M W 8.8 Maule earthquake in central Chile (Heki and Matsuo, 2010), and the 2011 M W 9.0 Tohoku earthquake in Japan (Matsuo and Heki, 2011).
The above studies on coseismic gravity change detection are focused on marine earthquakes, whereas studies of the far-field effect of gravity changes induced-by typical inland large earthquake are lacking. On April 25, 2015 at 06:11:26 UTC, a large, shallow earthquake (M W 7.8) occurred in the Gorkha region of central Nepal, rupturing a ~140 km segment of the central Himalayan front fault zone. Two aftershocks occurred on the same day at 06:15:22 (M W 6.7) and at 06:45:21 (M W 6.6), and on April 26 at 07:09:10 another aftershock (M W 6.7) occurred near the epicenter of the mainshock (USGS, 2015;Avouac et al., 2015). Since 2015, four SG stations (Lhasa, Lijiang, Wuhan and Beijing) have been operating in Chinese mainland, and the Wuhan and Lijiang stations are collocated with continuous high-frequency GPS observations. The aim of this paper is to detect far-field gravity changes caused by large earthquakes. These data might be useful for future studies in geophysics, such as improvement of dislocation theory. We investigated gravity changes related to the 2015 Nepal earthquake using SGs. First, we processed the raw gravity time series from SGs. Second, decimation filters were used to convert the raw data sampled at 1 s to 1 minute, and then the preliminary results of gravity changes were resolved. Third, the residual gravity signal at 1 minute was filtered to 3 hour intervals, consistent with the time resolution of near-surface hydrological models used in this study. Correction for this effect shows that significant far-field gravity changes were detected by the array of SGs at three stations (Lhasa, Lijiang and Wuhan station). Finally, we analyzed and interpreted the data in terms of VD and subsurface mass change by using cGPS and groundwater level (GWL) observation data, respectively.

SG Stations
The SG data used in this study were obtained from Lhasa, Lijiang, Wuhan and Beijing stations which are shown in Figure 1. The SG gravity station information is listed in Table 1. The OSG-057 at the Lhasa station, OSG-065 at the Wuhan station and OSG-066 at the Lijiang station are supported by Chinese Academy of Sciences; the iGrav-007 at the Wuhan station which was collocated with OSG-065 is supported by Wuhan University; and the iGrav-012 at the Beijing station is supported by National Institute of Metrology, China. SG records obtained at a 1 Hz sampling rate were used in this study.

Data Processing
Normally, the processing of SG data to obtain residual gravity, including a possible coseismic gravity change, is straightforward. Large accelerations due to high-frequency seismic waves tend to saturate the gravimeters, rendering the gravity data unusable at the time of the earthquake and shortly thereafter. Even if the gravity signal is not saturated, large horizontal accelerations can degrade the performance of the gravimeter (Imanishi et al., 2004). The recorded raw gravity data of the earthquake with 1 Hz sampling rate was deduced from theoretical Earth tides and ocean tide loading (Sun HP et al., 2019), air pressure, and polar motion effects using the program Tsoft (Van Camp and Vauterin, 2005). The results of residual gravity from OSG-065 on 25 April 2015 are shown in Figure 2. The residual gravity shows that the SG is partially out of range after the arrival of seismic waves during the mainshock, due to the measurement range of the SG covering ± 900 μGal. The results of the residual gravity from other SGs (from 22 April 2015 to 28 April 2015) were deduced similarly. Although each of the SGs was not saturated by seismic waves of the aftershocks, the DC levels of tilts in both the x and y directions did not appear to be biased by seismic waves of both mainshock and aftershocks, except for OSG-066 at the Lijiang station before the mainshock. In order to obtain an unbiased estimate of coseismic step, the span of the discarded data portion should be as short as possible.

Gravity Step Detection
To detect a coseismic gravity change, two periods of data were selected before and after the Nepal earthquake. Considering the earthquake released a large amount of seismic energy, the data were divided into two blocks. Block1 starts from 00:00:00 22 April to 06:10:00 25 April (about 10 minutes before seismic wave arrival of the mainshock), and to 02:30:00 24 April for OSG-066 at the Lijiang station due to the unbalanced tilt of the ball, which recovered again before the mainshock. Block2 starts from 22:00:00 25 April (about 16 hours after the mainshock) to 23:59:00 28 April. Continuous SG data with 1 Hz sample rate were corrected for body and oceanic tides, polar motion, and atmospheric effects as described above, then the residual gravity time series was lowpass filtered with a cutoff period of 0.00833 Hz and re-sampled to 1 min to suppress micro-seismic signals. The data from before and after the earthquake were separately fitted by a quadratic function to consider the different nonlinear trends: (1) Next, the gravity of the fitted functions at the last time of Block1 for the two blocks were calculated, for Block1 and for Block2. The gravity step can be determined by (Imanishi et al., 2004): Gravity steps were then computed, shown in Figure 3. The steps caused by the mainshock of the 2015 Nepal earthquake were very significant. Gravity is integrated, which means that a unique value includes the action of all the existing masses around the instrument, close by and very far away. Next, we will analyze the effects of VD and subsurface mass.

Vertical Deformation
The Wuhan and Lijiang SG stations are collocated with cGPS stations with 1 Hz and 30 s sampling rate, respectively. The cGPS data on 24, 25, and 26 April, 2015 were processed using the GIPSY software to remove phase and pseudorange data. The results are shown in Figure 4. The coseismic VDs were estimated by subtracting the averaged value on 24 April from 26 April (before and after the day of the Nepal earthquake). The results show that there are almost no permanent VD at both the Lijiang and the Wuhan stations. These results are consistent with the results provided by GNSS Data Products of the China Earthquake Administration (http://www.cgps.ac.cn). Considering the accuracy of cGPS and the VD of −3.9 ± 3.1 mm at the LHAZ station (Su XN et al., 2015) within several kilometers from OSG-057, the observed gravity step effect was about ~0.8 μGal, unaffected by the local or regional VD at the other three stations. However, the residual gravity signals were still affected by near-surface crust which consists of an un-saturated layer (soil-water layer), a saturated layer (groundwater, aquifers) and an impermeable layer.

Total Soil Moisture
Time series of gravity residual signals can be strongly influenced by hydrological processes, especially water mass variation. Therefore, hydrological effects should be modeled and removed so that the coseismic gravity changes can be appropriately estimated. The gravity variations due to hydrology can be separated into two major scales: local and regional (Llubes et al., 2004). The local scale is dominated by the Newtonian attraction of the underlying water masses. At the regional scale, surface and shallow water induces a global elastic deformation of the Earth, which has an effect on the gravity field through both mass redistribution and vertical movement.
The Global Land Data Assimilation System (GLDAS) NOAH V2.1 hydrological data with a 0.25° space resolution and a 3-hour time resolution (Rodell et al., 2004) can be used to simulate a variety of parameters, including soil moisture (at depths of 0−10 cm, 10− 40 cm, 40−100 cm and 100−200 cm), canopy water, and snow wa-  ter equivalent, among others. The products of GLDAS are freely available (https://disc.gsfc.nasa.gov/dataset?keywords=GLDAS). For the computation of the global hydrology effect, the regional effect was obtained by considering the parameters of loading and of Newtonian attraction at distances greater than 0.1° for convolution between the surface mass distribution and Green's functions. Due to the sensitivity of the ground gravimeters essentially limited to 1 km 2 around the instrument (Van Camp et al., 2017), water within 0.1° has a dominantly Newtonian effect. Thus the equivalent water height can be converted into gravity signal using the Bouguer conversion ratio of 0.42 μGal per mm of equivalent water height (Creutzfeldt et al., 2008). The equivalent water heights were extracted from GLDAS models by bilinear interpolation in terms of the four nearby grid values (Mikolaj et al., 2016).
The OSG-065 installed at Wuhan Observatory of Dynamic Geodesy is located at the top of Yanjia Mountain in the suburb of Wuhan. In addition to a GPS Trimble receiver, a groundwater level monitoring instrument in the borehole well and a meteorological instrument were also installed. The unconfined GWL is recorded with a ten-minute interval, with rainfall at a half hour interval. Although there is no effect on gravity resulting from the local vertical deformation before and after the 2015 Nepal earthquake, the combined co-site observations of the SG, cGPS, GWL and rainfall can be used comprehensively to monitor gravity due to large earthquakes by eliminating local, regional and global environmental perturbations. The solid line in Figure 5 shows the estimated gravity change associated with total soil moisture content from GL-DAS models in April, 2015. It can be seen that gravity increases steeply during precipitation events, and it deviates to a maximum of ~3 μGal according to the precipitation amount at each event. After precipitation, gravity decreased linearly. Gravity change shown inside the dashed rectangle indicates that there was no rainfall during this period and gravity decreased steeply with minor undulations. Therefore, the step of gravity change related to rainfall or the unsaturated layer at the Wuhan station can be excluded.
The results of total soil moisture re-distribution are plotted in Figure 6, showing that the gravity variations at each of the four stations ranges between −1.5 to 0.5 μGal, and the amplitudes are about 2.0, 1.0, 0.5 and 0.5 μGal at Wuhan, Lijiang, Lhasa and Beijing stations, respectively. There is no obvious step before and after the Nepal earthquake mainshock. To remain consistent with the time resolution of GLDAS models (3 hour intervals), the residual gravity time series with 1 min intervals observed by SGs at Lhasa, Lijiang and Wuhan station were downsampled to a 3 hour interval, with a cut-off period of 8 Circle Per Day (CPD). Correcting for the effects of toal soil moisture re-distributions, the residual gravity time changes at each station were derived, and the final gravity steps before and after the 2015 Nepal earthquake were separately fitted by a quadratic function. The estimated steps are shown in Figure 7 and Table 2.

Groundwater Level
The unconfined GWL at the Wuhan station with 10-min intervals is shown in Figure 8 Figure 5 at the four stations.

Earth and Planetary Physics
doi: 10.26464/epp2021018 145 ground surface level, and the GWLs were measured using pressure gauges at a depth of 10 m below the ground surface. After the 2015 Nepal earthquake, a step-like decrease in groundwater levels was detected. The observed total level change was −92.0 mm. Utilizing an effective porosity of 0.2 at the Wuhan station determined by Zhang WM et al. (2001), the gravity change due to local groundwater level change was about −0.8 μGal. Therefore, the gravity changes caused by the 2015 Nepal earthquake were 8.8 ± 0.2 μGal and 8.7 ± 0.1 μGal, detected by iGrav-007 and OSG-065, respectively.
Although there are no co-located groundwater wells at Lhasa, Lijiang and Beijing station, separate studies reported coseismic groundwater level changes of −2.0 mm and 16 mm at Lhasa geomagnetic station and Lijiang earthquake agency, respectively  Note: Contributions of toal soil moisture re-distributions were corrected for gravity steps of 3 hours. (Zhang B et al., 2015;Ren HW and Zhang L, 2017). Using an empirical porosity of 0.1, the gravity changes were ~0.1 μGal and about 0.1 μGal at Lhasa and Lijiang station, respectively, while the contribution of GWL changes can be ignored due to such a small amplitude.
In addition, there is a confined groundwater well with 1-day intervals at the Shisanling seismological station, less than 10 km from iGrav-012. The coseismic GWL change was 1.2 mm, obtained from the difference between April 24 and April 26 observations, yielding increased gravity values of ~0.1 μGal which can also be ignored.

Comparison with the Results of Dislocation Theory
We computed coseismic gravity changes using the spherical dislocation theory by Sun WK andOkubo (1993, 1998) with the USGS model (USGS, 2015). Results are depicted in Figure 9, which show that the general distribution exhibits a four-quadrant pattern, with positive change in the east-west direction and negative in the north-south direction. All four stations are located in the area of positive gravity changes but are less than 0.2 μGal, and even smaller at Lhasa and Beijing station which are located near the zero value of the contour. In general, the closer a station is to the epicenter, the larger the coseismic gravity change is. Although the Lhasa station is closer to the epicenter than Lijiang and Wuhan stations, the expected gravity change is smaller at the Lhasa station than the other two. Compared with the observed gravity changes at the four stations, three stations (Lijiang, Wuhan and Beijing) have the same signs as the theoretical changes, and one is opposite, but the observed magnitudes are significantly larger than the theoretical ones. Seismic dislocation theory plays a key role in seismic fault inversion and geodetic observation data interpretation. Dislocation theory is based on a regular geometric Earth model, such as semi-infinite space, homogeneous space, or layered sphere. The observed gravity changes from SGs do not coincide with the results of dislocation theory, which implies that a more realistic Earth model and reliable fault model are necessary (Zhang XL et al., 2016;Dong J et al., 2021).

Discussion and Conclusions
The main challenge of ground-based gravimetry is climate-induced mass changes such as precipitation, erosion, and changes of near-surface and subsurface water (Creutzfeldt et al., 2010;Van Camp et al., 2016). To detect the gravity step signal associated with a large earthquake, SG, with high-precision to ~0.1 μGal, provides a potentially powerful method to monitor mass changes in local and non-local regions. Therefore, hydrological, cGPS data and SG gravity data were collected in this study, and gravity changes were modeled using total soil moisture distributions from GLDAS models, VDs from cGPS, and GWL. The set scatter before and after the 2015 Nepal earthquake of OSG-056 at the Lhasa station was larger than the others; this may be due to its higher noise level (Zhang K et al., 2018) and shorter epicenter distance than the other SGs. Nonetheless, the gravity change caused by the earthquake can be properly fitted with an accuracy of 0.5 μGal.
Also, due to the lack of a co-site GPS at the Lhasa station, we could not estimate the precise result of coseismic VD; instead, results from other studies were used. The VD of −4.0 mm ranged from −7.0 mm to −0.8 mm, indicating a decreased trend which should account for ~0.8 μGal. Correction for this effect shows a significant ~−3.8±0.3 μGal decrease of the coseismic gravity change due to mass redistribution associated with the earthquake. Imanishi et al. (2009) pointed out that the coseismic gravity change observed by ground-based gravimetry consists of two effects: the change in Newtonian attraction due to the mass redistribution inside the Earth, and the VD of the Earth's surface. In general, these effects have comparable magnitude and show complicated dependence on the azimuth and the epicentral distance to gravity stations, and ground-based gravity observation alone cannot differentiate between the two effects.
The Nepal region also suffered three other aftershocks (~M W 7.0) besides the mainshock, and there was almost no obvious gravity change step observed by SGs using the same method as described above, so we can conclude that the influence of the aftershocks was minimal and can be ignored, since their magnitude is far smaller than the mainshock. The energy released indicates the destructive power of an earthquake. While the difference of magnitudes is 1.1 between M W 7.8 and M W 6.7 may seem minor, magnitude is a logarithmic scale, thus an increase of a unit of magnitude is about 32 times more energy (Richter, 1958), which means that the energy released by the Nepal mainshock is roughly 45 times that of the aftershock. . Distribution of coseismic gravity changes caused by the 2015 Nepal earthquake, as computed using the spherical dislocation theory (Sun WK andOkubo, 1993, 1998). Excluding the instrument offsets, especially the consistency of the two SGs (OSG-065 and iGrav-007) at the same place of the Wuhan station, reliable far-field gravity changes caused by the 2015 Nepal earthquake were found in this study, and the changes were −3.0 ± 0.5 μGal (OSG-057) at the Lhasa station, 7.3 ± 0.2 μGal (OSG-066) at the Lijiang station, 7.9 ± 0.1 μGal (OSG-065) and 8.0 ± 0.2 μGal (iGrav-007) at the Wuhan station, and 0.9 ± 0.2 μGal (iGrav-012) at the Beijing station.