A high-resolution time-variable terrestrial gravity field model of continental North China

Abstract

Given the spatial resolution of ~300 km in Gravity Recovery and Climate Experiment (GRACE) measurements, accurately quantifying mass variations at smaller scales proves challenging. Here, we present a high-resolution time-variable gravity field model of continental North China. This model, denoted as IGP-NorthChina2022TG, contains 15 gravity field solutions determined using an innovative approach that relies on terrestrial gravity measurements and Slepian basis functions. IGP-NorthChina2022TG provides degree 150 gravity changes (spatial resolution of ~120 km) on a semi-annual basis from September 2009 to September 2016, in contrast to the monthly degree 60 GRACE solutions. Despite their different temporal resolutions, the good agreement between GRACE and the ground-based results up to degree 60 confirms the robustness and reliability of the proposed method and favors the combination of these two types of measurements. The gravity changes with much finer spatial resolution from IGP-NorthChina2022TG could complement GRACE solutions for sub-regional scale investigations in North China.

Introduction

The time-variable gravity field, which reflects mass variations within the Earth system, is one of the keys to understanding the changing Earth. The mass variations are caused by the transport and redistribution of mass at the Earth’s surface and in its interior (crust, mantle, and core). Since its launch in 2002, the Gravity Recovery And Climate Experiment (GRACE) and its follow-on mission have revolutionized how we observe and understand global mass variations by measuring the time-variable gravity field¹.

North China is a region with well-developed active faults (Fig. 1a) and a high seismic hazard level. North China is also a substantial industrial base and grain-producing area, featuring a dense population and high water demand. Because of insufficient precipitation, groundwater has become the primary source of intensive industrial, agricultural, and domestic water use in North China and has been overexploited since the 1970s^2,3. Long-term groundwater overexploitation has caused the development of numerous cones of depression in the North China Plain (NCP; Fig. 1b), resulting in severe land subsidence (Fig. 1a) that affects sustainable economic and social development^4,5,6. Furthermore, overexploitation of groundwater in North China unloads the lithosphere and contributes to the deformation of bedrock, which may trigger the occurrence of earthquakes in this region^7,8. Time-variable gravity measurements from GRACE not only offer an independent tool to monitor and analyze groundwater storage (GWS) depletion in North China^{3,9,10,11,12,13,14,15,16,17}, but also facilitate estimates of the uplift rate (at the sub-millimeter per year level) in the same region^7,8,18.

Fig. 1: Tectonic setting and severe subsidence in North China.

a Green lines delineate the boundaries of the active blocks: Xiyu (A₁), Tibetan Plateau (A₂), South China (A₃), North China (A₄), and Northeastern Asia (A₅). White lines delineate the boundaries of the secondary active blocks: Ordos (A_4-1), North China Plain (A_4-2), and Eastern Shandong-Yellow Sea (A_4-3). All block data are taken from Zhang et al.⁷⁵. Thin red lines show the active faults⁷⁶. The GPS-derived vertical displacement rates in North China^8,17 are indicated by yellow and blue arrows, where negative rates represent subsidence. b Blue triangles show the locations of the repeated survey gravity stations. Magenta, green, black, and yellow dashed lines denote the boundaries of the region of interest, the North China Plain (NCP), the provincial borders, and a region with sparse stations, respectively. Number 1 denotes the Tianjin Municipality, and numbers 2 to 6 denote Shijiazhuang (capital city of Hebei Province), Taiyuan (capital city of Shanxi Province), Zhengzhou (capital city of Henan Province), Hefei (capital city of Anhui Province), and Jinan (capital city of Shandong Province), respectively. The shaded topography in the background is from SRTM15+⁷⁷.

Despite many successful applications in North China, accurately quantifying mass variations at small scales is still challenging due to the coarse spatial resolution (~300 km^3,19) of GRACE measurements. Moreover, for popular GRACE products in terms of spherical harmonics (SH), an additional postprocessing step such as decorrelation and spatial smoothing^20,21 is necessary to reduce the correlated and high degree errors. However, postprocessing is accompanied by signal attenuation and distortion^3,22,23. As alternatives to SH products, the mascon GRACE products provide grid results with improved resolutions that can be used directly without postprocessing^24,25,26,27. However, previous studies have pointed out that the spatial resolution of the mascon products is still relatively low, ranging from 200 km to 300 km^28,29.

Unlike the high-orbit GRACE observations, the terrestrial time-variable gravity measurements are collected at the Earth’s surface by regularly observing all stations in a network^30,31. For example, the field observations for the gravity network in North China (GNNC) are performed semiannually³². As part of the continental-scale gravity network established across the Chinese mainland, GNNC consists of 18 absolute-gravity stations and approximately 500 relative-gravity stations³². The distance between adjacent stations in GNNC ranges from several hundred meters to approximately 200 km. Consequently, gravity measurements from GNNC can provide information at a higher spatial resolution than the GRACE results. Moreover, terrestrial measurements are more accurate at short wavelengths due to the inherent signal attenuation with altitude of the gravity field³³.

In this study, we focus on developing a high-resolution time-variable gravity field model in North China solely from terrestrial measurements. The development includes 15 repeated survey campaign measurements from GNNC between September 2009 and September 2016 (see the Methods section, dataset subsection). Here, we present the derived time-variable gravity field model in North China with 15 solutions up to degree 150 (d150), as well as a tailored method for recovering regional gravity field based on Slepian basis functions^34,35,36 and terrestrial measurements (see Methods, Supplementary Methods 1 and 2). In addition, we show that the derived d150 solutions provide gravity changes with finer spatial resolution and stronger signal strength compared with GRACE-derived results for investigations at sub-regional scales.

Results and discussion

Overview of the modeling methods

We propose a tailored method for regional gravity field modeling using Slepian basis functions, which is comprehensively described in the Methods section, Supplementary Methods 1 and 2. This method consists of two key stages: model selection and the mitigation of broadband leakage. A brief overview of this method is outlined below. In the first stage, the maximum spherical expansion degree N and the optimal truncation level J_opt are determined using a two-step procedure. Within a given region, a subset of the Slepian basis functions tends to be well-concentrated, while the others are not. Consequently, it is reasonable to use the subset of well-concentrated Slepian basis functions to optimally model the signal confined to the specific region and omit the rest. The number of this subset of Slepian basis functions is termed the optimal truncation level and denoted by J_opt. Once these parameters (N and J_opt) are determined, their corresponding model is selected as the optimal one. In the second stage, the broadband leakage is mitigated using an iterative strategy based on Gaussian spatial averaging. This is an important issue in Slepian analyses that few studies have addressed. Broadband leakage arises from the discrepancy between the maximum expansion degree N and that contained in the input signals (which is infinity for in situ measurements). When modeling a degree N field, the high degree (>N) components in the input signal will affect the estimation of the coefficients at lower degrees (≤N). As detailed in the Methods section, Supplementary Methods 1 and 2, the proposed strategy effectively mitigates the broadband leakage.

Comparisons with GRACE-derived gravity changes

Before moving on to the d150 gravity field recovery in North China, we validated the proposed method (see Methods) by recovering a series of degree 60 (d60) fields and comparing them with the GRACE-derived results. Using terrestrial gravity measurements, we first set the optimal truncation level J_opt and the Gaussian spatial averaging radius r [km] for the d60 field recovery to 3 and 512, respectively. A series of 15 d60 gravity solutions were obtained using these parameters and the terrestrial gravity measurements (see Supplementary Note 1 and Supplementary Figs. 7 and 8). We then estimated the rates of the d60 gravity changes ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ (where TG indicates terrestrial gravity) in North China by least-squares fitting a linear model to the time series at each grid node of the 15 d60 solutions. The seasonal variations are not estimated due to the semi-annual temporal resolution of TG measurements (see Datasets). Figure 2a shows the derived linear rates and reveals an overall increasing pattern over the entire region of interest, with a dominant gravity increase centering around the city of Shijiazhuang in the northern part of the region of interest. Next, we examined the vertical displacements in the region of interest. The GRACE results are derived at the space-fixed points^37,38,39, while the terrestrial measurements are derived on the Earth’s surface, which is associated with vertical displacements (either uplift or subsidence). For example, the NCP (see Fig. 1b and Supplementary Fig. 9b) in the region of interest has suffered from severe land subsidence during the past several decades^5,6,13,17. Therefore, we need to remove the gravity effects induced by the vertical displacements in the region of interest prior to the comparison. It is accomplished by incorporating the GPS vertical displacement measurements (Supplementary Note 2). We recovered a d60 field, ${{{{{\rm{d}}}}}}{h}_{{{{{{\rm{G}}}}}}{{{{{\rm{PS}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ (Supplementary Fig. 9c), from these GPS vertical displacement rates using the same above-determined J_opt and r parameters.

Fig. 2: Comparison of terrestrial- and GRACE-derived gravity rates in North China.

a Rates of the degree 60 gravity changes derived using Slepian functions and terrestrial gravity measurements; b rates of the gravity changes induced by vertical displacement up to degree 60; c residual rates (subtracting panel b from panel a); d ensemble rates of Slepian-derived GRACE gravity changes; e ensemble rates of DDK5 (denoising and decorrelation kernel version 5)-derived GRACE gravity changes; and f ensemble rates of Gaussian smoothing-derived GRACE gravity changes. The magenta line denotes the boundary of the region of interest. The green line shows the boundary of the NCP, the black line indicates the land-ocean boundary and the gray lines show the provincial borders.

We then accounted for the gravity effects induced by vertical displacements using different gravity gradient ${\gamma }_{g}$ in North China (see Methods), namely the free-air gradient (${\gamma }_{g}$ = –3.086 μGal cm^–1), the Bouguer gradient (${\gamma }_{g}$ = –1.967 μGal cm^–1) and the Bouguer-corrected free-air gradient (BCFAG_spherical) of a spherical source (${\gamma }_{g}$ = –3.086 + 4${{{{{\rm{\pi }}}}}}G{\rho }_{0}$/3 μGal cm^–1)^40,41, where G is the gravitational constant and ${\rho }_{0}$ is the density of the source. The BCFAG_spherical represents the gravity changes associated with a change in elevation, taking into account the mass accompanying the deformation³³. Assuming a spherical magma source, the BCFAG_spherical is widely adopted in dealing with gravity changes in volcanic systems^33,40,41. Since the land subsidence in the NCP has a strong relationship with groundwater overexploitation⁵, we assume a spherical source of groundwater by analogy with the volcanic case when applying the BCFAG_spherical correction. With groundwater density of 1000 kg m^–3, the BCFAG_spherical results in ${\gamma }_{g}$ = –2.806 μGal cm^–1.

The gravity effects induced by vertical displacements ${\gamma }_{g}{{{{{\rm{d}}}}}}{h}_{{{{{{\rm{GPS}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ with ${\gamma }_{g}$ as the free-air and Bouguer gradients (see Eq. 4 in Methods) are shown in Supplementary Figs. 12a and 12d, respectively. The gravity effects ${\gamma }_{g}{{{{{\rm{d}}}}}}{h}_{{{{{{\rm{GPS}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ with ${\gamma }_{g}$ as the BCFAG_spherical gradient are shown in Fig. 2b. After accounting for the vertical displacement effects, the residual gravity rate ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ (see Eq. 4) in Fig. 2c is comparable to the GRACE-derived gravity rate.

As an independent estimation of the gravity changes, three sets of GRACE-derived gravity rates (see Supplementary Note 3 for details) over the region of interest are presented in Fig. 2d–f. These d60 rates using GRACE data all show a consistent pattern of gravity decreases over the region of interest, with the most notable decrease close to the border of three provinces (Hebei, Henan, and Shandong). The rates with the least attenuated signal strength among the GRACE-derived d60 rates are those derived using Slepian functions (Fig. 2d), which do not require smoothing or filtering. The d60 results obtained using the DDK (Fig. 2e) and Gaussian smoothing filters (Fig. 2f) further attenuate the signal power. For example, the mean d60 rates in the NCP from 2009 to 2016 are −1.27 (Slepian-derived GRACE, Fig. 2d), −1.00 (DDK-derived GRACE, Fig. 2e), and −0.99 μGal yr⁻¹ (Gaussian-derived GRACE, Fig. 2f). Hence, we use the Slepian-derived GRACE results in Fig. 2d as a benchmark in comparison with the d60 rates derived using terrestrial measurements. Note that parts of the NCP are outside the region of interest (Fig. 1b). Hereafter, we refer to the majority of the NCP, which lies inside the region of interest, as NCP for brevity.

Using different gravity gradients (Eq. 4), we derived three sets of d60 gravity residual rates ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ based on the terrestrial gravity and GPS measurements in North China, as shown in Supplementary Figs. 12b, 12e and 12h (Supplementary Note 4), respectively. The differences between terrestrial- and GRACE-derived gravity rates (Fig. 2d) in North China are shown in Fig. 3 and Supplementary Figs. 12c, 12f and 12i. In the region of interest, the correlation coefficient between ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ derived using free-air gradient and GRACE rates has the largest value (0.898), followed by that between ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ using the BCFAG_spherical and GRACE rates (0.878), and the correlation coefficient between ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$ derived using the Bouguer gradient and GRACE rates has the smallest value (0.671). These correlation coefficients and the differences in Fig. 3 indicate that the Bouguer gradient is inadequate to account for the gravity effects induced by vertical displacements in North China, especially in the NCP. Over the region of interest, the difference in the free-air case (Fig. 3a) and that in the BCFAG_spherical case (Fig. 3b) show negative-positive-negative spatial patterns from north to south. In comparison, the difference in the Bouguer case (Fig. 3c) shows a prominent positive pattern. Moreover, we examined the mean and root mean square for the differences of the three cases in the NCP. The mean and root mean square are –0.040 and 0.147 μGal yr^–1 for the BCFAG_spherical case, –0.365 and 0.379 μGal yr^–1 for the free-air case, and 0.909 and 0.962 μGal yr^–1 for Bouguer case. These analyses show that the BCFAG_spherical performs the best in the region of interest among the three vertical gravity gradients.

Fig. 3: Differences between terrestrial- and GRACE-derived gravity rates in North China.

These differences are calculated by subtracting GRACE-derived rates in Fig. 2d from terrestrial residual rates corrected by a Bouguer-corrected free-air gradient of a spherical source; b free-air gradient; and c Bouguer gradient. The magenta lines denote the boundary of the region of interest, the green lines show the boundary of the NCP, the black lines indicate the land-ocean boundary, and the gray lines show the provincial borders. Note that the color bar is different from that in Fig. 2.

This good agreement between the satellite- and ground-based results confirms the robustness and reliability of the model selection and leakage mitigation strategies (see Methods). Moreover, this comparison implies that ground-based gravity field solutions can be applied to validate satellite-based solutions. In addition, the signals derived using Slepian functions (Fig. 2c and d) are well concentrated in the region of interest. Those derived using the DDK and Gaussian smoothing filters are smeared out, inducing spatial leakage (evident at the land-ocean boundary). This contrast further highlights the benefits of Slepian functions for regional analyses.

Degree 150 gravity field solutions

Based on the proposed methods, we recovered a series of 15 solutions up to degree 150 (d150) solely from repeated campaign measurements at the 287 gravity stations, covering the period from September 2009 to September 2016 (Fig. 4a–o). The spatial resolution is ~120 km (half-wavelength). We call these solutions IGP-NorthChina2022TG, where TG indicates that each solution was derived solely from terrestrial gravity measurements.

Fig. 4: The derived time-variable gravity field model (degree 150) in North China.

Spatial pattern of the degree 150 solutions in a September 2009; b April 2010; c September 2010; d April 2011; e September 2011; f April 2012; g September 2012; h April 2013; i September 2013; j April 2014; k September 2014; l April 2015; m September 2015; n April 2016; and o September 2016. The mean of all 15 solutions is removed from each solution. This model is called IGP-NorthChina2022TG (where TG indicates terrestrial gravity) and has a temporal coverage from September 2009 to September 2016 on a semi-annual basis. The magenta line denotes the boundary of the region of interest, the black line indicates the land-ocean boundary, and the gray lines show the provincial borders.

Next, we estimated the linear rates from September 2009 to September 2016 for the degree 150 gravity changes, ${g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$, in North China from IGP-NorthChina2022TG (Fig. 4) using least-squares fitting, as we did for ${g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}60}/{{{{{\rm{d}}}}}}t$. The derived rates, ${g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ (Fig. 5a), show much finer spatial resolution and stronger signal strength compared with the rates from GRACE, revealing a decreasing–increasing–decreasing pattern from north to south over the region of interest. The first decreasing pattern is to the northwest of Beijing, with a negative peak reaching approximately −4 μGal yr⁻¹. Then, there is an increasing pattern over the middle part of the region of interest, including the majority of Shanxi Province, the southern part of Hebei Province, the western part of Shandong Province, and the northern part of Henan Province (see Fig. 1b for the locations of these provinces). The region with the most notable gravity increase (approximately 10 μGal yr⁻¹) is located around Hengshui, close to the border of the Hebei and Shandong provinces. This finding is similar to previous terrestrial gravity measurement studies, e.g., results from Shen et al.⁴² from 2009 to 2013 and Hao et al.⁴³ from 2010 to 2014. As part of the NCP, this region also suffers from groundwater depletion and severe land subsidence^6,8,44,45. Finally, there is an elongated, northwest–southeast-oriented, decreasing pattern over the southern part of the region of interest with two peaks (approximately −5 μGal yr⁻¹). One peak is just south of Zhengzhou in Henan Province, while the other is near Hefei in Anhui Province.

Fig. 5: Comparison of d150 gravity rates.

a Gravity rates from the IGP-NorthChina2022TG degree 150 model (September 2009–September 2016), yellow crosses denoting the locations of the absolute gravity (AG) stations; b gravity rates induced by vertical displacement up to degree 150; c rates of hydrological effects from the WaterGAP Global Hydrology Model (WGHM) up to degree 150; d rates of the degree 150 residual gravity changes derived by subtracting panel b from panel a; e rates of the degree 150 residual gravity changes derived by subtracting panel c from panel d; and f difference between panel d and those in Fig. 2d. The yellow-red dashed line denotes the boundaries of the active blocks⁷⁸, with S₁ denoting the Taihangshan sub-block, S₂ denoting the Hebei-Shandong sub-block, and S₃ denoting the Henan-Huai sub-block. S₁, S₂, and S₃ belong to the NCP block (A_4-2 in Fig. 1a). The magenta lines denote the boundary of the region of interest, the green lines show the boundary of the NCP, the black lines indicate the land-ocean boundary, and the gray lines show the provincial borders. Note that the color bar is different in panel c.

Both terrestrial and GRACE measurements capture integrated gravity changes originating from geophysical phenomena occurring at the surface and inside the Earth. Hence, the degree 150 gravity rate, ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$, is also an integration representing the mass changes from those phenomena. To help identify the gravity changes beneath the Earth’s surface, we further removed gravity rates induced by vertical displacements and hydrological signals from the integrated rate, ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$. We first converted the d150 GPS vertical displacement rate, ${{{{{\rm{d}}}}}}{h}_{{{{{{\rm{G}}}}}}{{{{{\rm{PS}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ (Supplementary Fig. 9d), into gravity changes using the BCFAG_spherical gradient (see Methods); the result is shown in Fig. 5b. Then, we used output from the WaterGAP Global Hydrology Model (WGHM) 2.2d^46,47 for the same months as the terrestrial gravity measurements to estimate the gravity effects induced by hydrological signals. The WGHM output includes changes in canopy storage, soil moisture storage, water-equivalent snow storage, and surface water storage (i.e., lakes, reservoirs, rivers, and wetlands). The sum of these storage changes in each month is first converted to gravity changes, from which we then estimate the rates of the hydrological signals, ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{Hydro}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$. As shown in Fig. 5c, the variations in the derived ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{Hydro}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ from WGHM 2.2d are much smaller than those in the gravity rate, ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$. For example, the mean ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{Hydro}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ (−0.05 μGal yr⁻¹) only accounts for ~5% of the mean ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ over the region of interest (0.93 μGal yr⁻¹). We estimated two sets of degree 150 residual rates, ${g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ (see Eq. 4 in Methods), from September 2009 to September 2016 in North China. The first set (Fig. 5d) was derived by removing the linear rates induced by vertical displacements from the gravity rates ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ (Fig. 5a). The BCFAG_spherical gradient was adopted in this procedure. The second set (Fig. 5e) was derived by further removing the hydrological rates ${{{{{\rm{d}}}}}}{g}_{{{{{{\rm{Hydro}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ from the results in Fig. 5d. The first set of ${g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ could be directly compared with the GRACE-derived rates in Fig. 2d, since both results are based on observations and not corrected for hydrological contributions. The second set of ${g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ aims to isolate the effects primarily caused by GWS changes and tectonics. Note that any uncertainty of the adopted model for estimating the hydrological contributions will contribute to the uncertainty of the d150 residual rates. Here, the second set of ${g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$ is affected by the uncertainty and limitation of the WGHM-based water storage changes in canopy, soil, snow, lakes, reservoirs, rivers, and wetlands. As we show later, WGHM 2.2d does not include all large reservoirs in North China in its estimation of surface water storage (Supplementary Fig. 14c). Previous study also shows that the WGHM-derived interannual rates are with large uncertainty⁴⁸.

Figure 5d, e show similar spatial patterns and roughly reveal a decreasing–increasing–decreasing pattern from north to south over the region of interest. A decreasing pattern is concentrated over the NCP with two peaks, with one peak (approximately −6 μGal yr⁻¹) around Beijing and the other (approximately −4 μGal yr⁻¹) around Hengshui. This pattern highlights a sub-regional difference in the decrease in the NCP, indicating that the decrease is faster in the eastern regions than in the piedmont regions. An increasing pattern occurs over the majority of Shanxi Province, with its peak (approximately 4 μGal yr⁻¹) close to the border between Shanxi and Henan provinces. This positive rate is clearly separated from the negative rate over the NCP along the Taihang Mountains (see Fig. 1b). A decreasing pattern with northwest-southeast orientation occurs over the southern part of the region of interest afterward. This pattern is similar to the negative pattern of the d150 gravity rates, ${g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{d}}}}}}150}/{{{{{\rm{d}}}}}}t$, in the same region but with a slightly stronger signal strength (negative peak value of approximately −6 μGal yr⁻¹) and a negative peak in Anhui Province moving northward.

We then compared the first set of d150 residual rates (Fig. 5d) with the d60 rates from GRACE (Fig. 2d). First, the d150 rates in Fig. 5d show a larger amplitude. For instance, the mean gravity rate in the North China Plain in Fig. 5d is –2.18 μGal yr⁻¹ (with a root mean square of 2.74 μGal yr⁻¹), nearly two times as large as the estimate of GRACE (–1.27 μGal yr⁻¹ in Fig. 2d). Second, the d150 rates reveal more spatial details in North China. In principle, the d60 GRACE rates capture an overall decreasing trend in the region of interest. After removing these d60 rates from the first set of d150 rates, we obtained the degree 61 to 150 spherical harmonic components of gravity rates from September 2009 to September 2016 in North China. Prominent variations are visible to the northwest of Beijing (decreasing pattern), the south of Shanxi Province (increasing pattern) and the south of Henan Province and the northwest of Anhui Province (decreasing pattern) in the region of interest in Fig. 5f; however, these variations are not captured by GRACE.

Relationship between the GWS-induced gravity rates and the estimated d150 residual rates

We collected GWS statistics from the annual Water Resources Bulletins (WRB) from two municipalities (Beijing and Tianjin) and five provinces (Hebei, Henan, Shandong, Shanxi, and Anhui). For the WRB of each province, the GWS statistics are given city by city. All the GWS data provided in the WRB are volumes (units in 10⁸ m³). The GWS statistics for Beijing, Tianjin, Hebei, Henan, and Shandong cover the years from 2009 to 2016, the same as our study period. However, the available GWS statistics for Shanxi and Anhui do not completely cover the study period, with data in Shanxi from 2011 to 2016 and data in Anhui from 2009 to 2013 (with missing data in 2010). Figure 6a shows the rates of the GWS volume changes estimated from these official statistics for each city in the region of interest. After dividing the area for each city, we converted the rates of the volume changes into rates of equivalent water height changes; the result is shown in Fig. 6b. These water height change rates are then converted to d150 gravity changes using Slepian functions (see Supplementary Methods 1 and 2 for details), as shown in Fig. 6c.

Fig. 6: Water Resources Bulletin– and WGHM-derived d150 gravity rates.

a Spatial pattern of shallow groundwater storage (GWS) rates in North China based on the official water resources bulletins from two municipalities (Beijing and Tianjin) and five provinces (Hebei, Henan, Shandong, Shanxi, and Anhui); b rates of the equivalent water heights (EWHs) inferred from panel a; c d150 rates of gravity changes from the equivalent water height changes; and d d150 GWS rates from WGHM. The magenta lines denote the boundary of the region of interest, the green lines show the boundary of the NCP, the black lines indicate the land-ocean boundary, and the gray lines show the provincial borders. Note that each panel has an individual color bar.

The WRB-based equivalent water height and gravity rates (Fig. 6b, c) highlight the depletion of GWS in the NCP and Henan Province, as well as the increase in GWS in Shanxi Province and part of Henan Province. It is worthwhile to note that possible GWS changes in deep aquifers are not included in Fig. 6b, c.

The WRB-based d150 rates indicate that (1) the most apparent depletion of GWS in the NCP occurred in the eastern regions from 2009 to 2016, while the piedmont regions in the NCP experienced a GWS increase during the same period; (2) two NCP regions show similar depletion rates, one occurring around Beijing and the other occurring to the south of Hengshui. The majority of Henan Province also experienced GWS depletion, with the maximum depletion occurring in Luoyang and Pingdingshan; (3) GWS increases occurred throughout the majority of Shanxi Province, with the maximum increase covering the cities of Jinzhong, Changzhi, and Jincheng. GWS increases with comparable amplitudes are found in Zhumadian and Xinyang, in the southeastern part of Henan Province. Except for the increase in Henan Province, the spatial pattern of the d150 WRB-based GWS rates (Fig. 6c) is in good agreement with the d150 ground-based residual rates (Fig. 5e).

Despite their similar spatial patterns, there are large differences in magnitude between the WRB-based GWS rates (Fig. 6c) and the ground-based residual rates (Fig. 5e). For example, the mean d150 rate in the NCP induced by groundwater depletion based on the WRB is −0.07 μGal yr⁻¹ from 2009 to 2016. Conversely, the mean d150 ground-based residual rate in the NCP is −2.07 μGal yr⁻¹ from September 2009 to September 2016. The WRB-based GWS rate accounts for less than 3% of the ground-based residual rate. This discrepancy may be explained by the fact that (1) their temporal coverages are different. The official bulletins are issued annually; therefore, the WRB-based rates cover a longer period (January 2009–December 2016). Recent studies have revealed a sub-regional recovery of GWS in the NCP after 2015, resulting from the South-to-North Water Diversion Project^15,16. This recovery has alleviated the GWS depletion in the NCP, resulting in a lower GWS rate from January 2009 to December 2016 than from September 2009 to September 2016; (2) the WRB-based estimate only includes GWS changes in the shallow unconfined aquifers^3,11. In contrast, the ground-based residuals include GWS changes from both shallow unconfined and deep confined aquifers, as well as gravity changes induced by possible tectonic effects. For example, the ground-based residual rate map highlights a region with an evident gravity increase in Shanxi Province (Fig. 5e). This region is on the Taihangshan sub-block (part of the North China Plain block), where possible upwelling of asthenospheric material has been reported⁴⁹. However, determining the effect of tectonic activities on gravity changes is beyond the scope of the current study; (3) the d150 ground-based residual rate may be subject to large uncertainty. The isolation of GWS depends on the removal of water storages using a hydrological model (WGHM in our case), whose uncertainty is unknown and will propagate into the uncertainty of the residual estimate.

The GWS estimates from WGHM 2.2d in North China were previously reported to be not sensitive to the actual GWS changes⁵⁰. Here, we estimate the d150 gravity rates from the GWS changes provided by WGHM 2.2d. The WGHM GWS changes for the same months as the terrestrial gravity measurements were used. The WGHM-based d150 GWS rates (Fig. 6d) reveal gravity decreases mainly confined to the NCP and the northern part of Henan Province. The mean rate in the NCP is −1.57 μGal yr⁻¹, smaller than that of the d150 ground-based residuals (Fig. 5e). The spatial patterns of the WGHM-based GWS rates differ greatly from those of the d150 residual rates: (1) the decrease in the NCP is more severe in the piedmont regions than in the eastern regions and (2) the center of the decrease in Henan Province is moving northward, and no decrease is found in Anhui Province. Although the d150 ground-based residual rate in Fig. 5e may be subject to large uncertainty, these differences imply that WGHM 2.2d may not fully capture the complete water storage changes in North China.

Acceleration in the gravity increase over the NCP

The mean d150 gravity changes averaged over the NCP show an overall increasing trend during the study period (Supplementary Fig. 13) and highlight an acceleration in the gravity increase from 2012 to 2016. This acceleration started in April 2012 with a rate of 3.52 ± 2.81 μGal yr^–1 from April 2012 to September 2016, nearly tripling the rate prior to April 2012. From April 2012 through April 2015, the acceleration rate reached 12.76 ± 3.24 μGal yr⁻¹. Zhao et al.¹⁷ reported an acceleration in land subsidence in the NCP during the period of 2013–2016 based on GPS vertical measurements, ranging from ~2 cm yr⁻¹ to ~10 cm yr⁻¹. Accordingly, the acceleration in the gravity changes from 2012 to 2016 is likely linked to the accelerated land subsidence. Note that the gravity changes show a decreasing trend from 2015 to 2016. If this decreasing trend continues, the acceleration in the gravity changes will be mitigated.

Gravity increases in Shanxi Province

An increasing pattern is evident in the d150 residual gravity rates in Shanxi Province (Fig. 5d). This pattern was also detected by Feng et al.⁵⁰ using the high-resolution Tongji-RegGrace2019 model, which is a regularized GRACE solution up to degree 180⁵¹. The d150 residual rates increase from north to south in Shanxi, and the region with the largest gravity increase (approximately 3.6 μGal yr⁻¹) covers the majority of the cities of Changzhi and Jincheng and part of the city of Linfen. The mean d150 time series of the gravity changes within this region (Supplementary Fig. 14a) shows an overall increasing trend, with a large drop occurring in 2012 and a dramatic increase from April 2013 to April 2015. The annual precipitations in Changzhi, Jincheng, and Linfen, according to the Shanxi WRB, show no dramatic changes from 2013 to 2015 (Supplementary Fig. 14b), implying that the gravity increase here is unlikely caused by surplus rainfall. The GWS rates derived from the WRBs are also small (approximately 0.25 μGal yr⁻¹; Fig. 6c) in this region. The region with an evident gravity increase also includes a large reservoir. The Zhangfeng Reservoir is located on the Qinhe River, a tributary of the Yellow River. The total capacity of the Zhangfeng Reservoir is 0.394 km³, approximately 100 times smaller than that of the Three Gorges Reservoir (~39.3 km³). The Zhangfeng Reservoir began impoundment in late 2007. The water level reached its maximum operational height of 759 m in October 2017, creating a lake with an area of 14.36 km² (http://www.chinawater.com.cn/newscenter/df/sx/201710/t20171020_493200.html). Even though WGHM provides estimations of water storage changes in reservoirs, those for the Zhangfeng Reservoir are not included in the latest WGHM 2.2d (Supplementary Fig. 14c).

Based on the published water level and volume records for the Zhangfeng Reservoir, we derived a quadratic relationship between the reservoir volume and water level (see Supplementary Note 5 for details). Applying a least-squares fitting to the water levels in Table S1, we obtained a water level change rate of 2.86 m yr⁻¹ at the Zhangfeng Reservoir from 2009 to 2016. Assuming that the area of the Zhangfeng Reservoir remains the same (14.36 km²), the mass changes are only associated with the changes in the water level. Note that this assumption leads to overestimating the gravity changes induced by water impoundment. Following this assumption, we converted the water level change rates into gravity rates (Supplementary Fig. 14f). The estimated gravity rates induced by the water impoundment of the Zhangfeng Reservoir are less than 0.10 μGal yr⁻¹, too small to be a substantial contributor to the gravity increase in the southern part of Shanxi Province.

Conclusions

This study aims to provide new information to the scientific community, i.e., the series of high-resolution time-variable gravity field solutions recovered solely from terrestrial gravity measurements (IGP-NorthChina2022TG). One key challenge arises when recovering the gravity field in North China using terrestrial measurements: the recovery of the regional gravity field from the sparsely (and unevenly) distributed and locally available observations. Terrestrial measurements are subject to geographical and environmental constraints, resulting in an uneven distribution of observations. For example, the distribution of observations is dense in the Beijing-Tianjin region, but sparse in Henan Province (Fig. 1b). In this study, we relied on a Slepian localization analysis to recover the regional gravity field in North China from terrestrial gravity measurements. We presented a tailored method for selecting the optimal parameters in the Slepian analysis. The good agreement between the GRACE- and ground-based results up to degree 60 confirmed the robustness and reliability of the proposed method, favoring the combination of these two types of measurements in the future. IGP-NorthChina2022TG also agrees reasonably well with independent in situ absolute gravity measurements.

IGP-NorthChina2022TG provides degree 150 (spatial resolution of ~120 km) time-variable gravity field solutions from September 2009 to September 2016 on a semi-annual basis. The analysis of the recovery process may provide useful information for refining the gravity network in North China (see Supplementary Note 6 for details). The degree 150 residual gravity rates estimated from IGP-NorthChina2022TG, with the gravity rates induced by vertical displacements removed, show much finer spatial resolution and stronger signal strength compared with GRACE-derived rates for investigations at sub-regional scales. These residual rates may complement GRACE-based analyses of GWS changes and other geophysical processes in North China. The degree 150 gravity time series also provide insights into the time-variable changes, e.g., the acceleration of the gravity increase after April 2012 in the NCP corresponds well with the acceleration of the land subsidence. As the gravity campaign over the network in North China continues and adjustment techniques evolve^31,52, regular refinement and updates of IGP-NorthChina2022TG will be possible.

Methods

Datasets

The core datasets used in this work are 15 repeated survey campaign measurements from a gravity network in North China between September 2009 and September 2016, provided by the Gravity Network Centre of China. Using multiple gravimeters such as Scintrex CG-5, LaCoste-Romberg G, and FG-5, these measurements were carried out roughly on a semi-annual basis, with one campaign around April and the other one around September in each year³². After correcting for effects from the earth tide, ocean tide, air pressure, and polar motion, the network adjustment was performed using the pyBACGS program^31,52, a part of the Geoist package (https://github.com/igp-gravity/geoist). The Bayesian adjustment model in pyBACGS allows for an effective estimation of the drift rate and the scale factor for each gravimeter, especially in the case of long-term gravity survey campaigns with multiple gravimeters⁵³. The network adjustments for all 15 campaigns yielded an average accuracy of 7.2–9.0 µGal. The locations of the 287 gravity stations covered by all 15 survey campaigns are shown in Fig. 1b. The positions of the most marginal stations generate the boundary (the magenta line in Fig. 1b) of the region of interest.

The auxiliary datasets include the Gravity Recovery and Climate Experiment (GRACE) gravity field solutions and the WaterGAP Global Hydrology Model (WGHM). We collected GRACE level-2 release 06 products from three processing centers: the Center for Space Research at the University of Texas, the GeoForschungsZentrum, and the Jet Propulsion Laboratory. Each level-2 product consists of a set of spherical harmonic (SH) coefficients up to degree and order 60. WGHM simulates water flows and water storage at a global scale. The spatial resolution of the WGHM grid, λ_WGHM, is 0.5°. The climate forcing for WGHM includes air temperature, precipitation, downward shortwave radiation, and downward longwave radiation, which are taken from a number of global meteorological datasets⁴⁷. Note that WGHM provides an estimate of the groundwater storage (GWS), as well as storage for nine other compartments (canopy, soil moisture, snow, local lake, global lake, local wetland, global wetland, reservoir, and river).

Model selection

The Slepian approach requires a region Γ, a maximum spherical harmonic expansion degree N, and the number J_opt of the Slepian basis functions (Supplementary Methods 2, Supplementary Figs. 1 and 2). The region Γ is the region of interest, which is outlined in Fig. 1b. The limiting factor for the maximum expansion degree N comes from the spatial data gaps⁵⁴, which, in our case, reflect the sparse spatial distribution of the gravity measurements. For the truncation level J_opt, even though the Shannon number provides an approximate estimate, an iterative selection procedure is often adopted in practical applications^54,55,56,57. Following the strategies devised by Plattner and Simons⁵⁴ and Han et al.⁵⁷, we propose a two-step procedure to select the optimal truncation level J_opt. The first step covers the strategy to find the number J_opt for a given expansion degree n. In the second step, we compare results for a series of expansion degrees n and select the maximum degree N. The strategy of step one is as follows.

(1) Construct synthetic gravity fields based on the model outputs from WGHM 2.2d. First, we retrieved the total water storage (TWS) by summing all the storage from the 10 compartments. We obtained 15 TWS results between September 2009 and September 2016, consistent with the period of the terrestrial gravity measurements. We then converted the TWS into equivalent water heights (EWHs). Next, we expanded these EWHs into SH coefficients ${f}_{{nm}}$. The SH coefficients ${f}_{{nm}}$ are up to degree and order 360 (the value 360 equals 180/λ_WGHM). Finally, we computed the gravity changes from these SH coefficients as the synthetic gravity fields based on the following relation:⁵⁸

$$\Delta g\left(\theta ,\lambda \right)=\frac{{GM}}{{a}^{2}}\mathop{\sum }\limits_{n=1}^{{N}_{0}}\left(n-1\right)\mathop{\sum }\limits_{m=-n}^{n}{f}_{{nm}}{Y}_{{nm}}(\theta ,\lambda )$$

(1)

where ($\theta ,\lambda$) denotes the colatitude and longitude of the point ($0\le \theta \le {{{{{\rm{\pi }}}}}}$ and $0\le \lambda < 2{{{{{\rm{\pi }}}}}}$), ${Y}_{{nm}}\left(\theta ,\lambda \right)$ are the spherical harmonic function of degree n and order m, GM is the product of the gravitational constant and the mass of the Earth, a is the semi-major axis of the Earth, and N₀ denotes different SH expansion degrees (N₀ = 360, 180, 150, 120, 90, and 60). The unit of the derived gravity changes is µGal (1 µGal = 1 × 10⁻⁸ m/s²).

(2) Generate true gravity signals for a given degree N₀ at the repeated survey gravity stations (see Fig. 1b) in the region of interest. After obtaining the synthetic gravity fields from September 2009 to September 2016 for different expansion degrees N₀ in step (1), we removed the mean from each N₀ series. All the tests were conducted using the gravity changes for different N₀ in September 2009. The derived gravity changes for different degrees N₀ in September 2009 are shown in Supplementary Figs. 3b–g. For a given degree N₀, the values at the repeated survey gravity stations linearly interpolated from the corresponding N₀ field are used as the true gravity signals in the following analysis.

(3) Perform the perturbation process for five alternative expansion degrees N₀ (180, 150, 120, 90, and 60). For each degree N₀, we first generated 100 realizations of Gaussian random noise with zero mean and a standard deviation of 3 µGal, the typical noise level for the relative gravity measurements⁵⁹. By adding these noise fields to each true N₀ gravity signal, we created 100 new datasets with perturbations ready for gravity field recovery using the Slepian approach. We recovered the gravity field for each new dataset using the J = 1, 2, 3,…, 50 Slepian basis functions. In total, we completed 5000 recoveries per given spherical harmonic degree N₀. The computation load is extensive, requiring approximately 12 h and 32 h to finish 5000 recoveries for N₀ = 150 and N₀ = 180, respectively, on a modern desktop computer (Intel i9 2.5 GHz, eight cores processor). Next, we adopted the Bayesian information criterion (BIC) that includes a penalty term increasing with the number of parameters in the model selection process. Given a set of Slepian-recovered model results with different truncation levels J, we have

$${{{{{\rm{BIC}}}}}}=J\cdot {{{{\mathrm{ln}}}}}\left(P\right)-2\,{{ln}}\left(L\right)$$

(2)

where P is the number of observations, and L is the maximum value of the likelihood function. The BICs were first evaluated over the entire region of interest (ROI) for different truncation numbers J. This series of BICs is denoted BIC_ROI. To highlight the influence of the sparsity of the gravity measurements, we also evaluated the BICs over the southern part of the region of interest, where the gravity stations are sparse (the region outlined by the yellow dashed line in Fig. 1b), resulting in another series of BICs denoted BIC_sparse. For each new dataset, we obtained two series of BICs (J = 1, 2, 3,…,50). Finally, we averaged the BIC_ROI and BIC_sparse series from all 100 datasets to obtain two mean BIC series. In these two derived mean BIC series, the truncation number J corresponding to the lowest BIC value is J_opt. The mean BIC_ROI and BIC_sparse series for N₀ = 150 are shown in Supplementary Fig. 4a. Both series first decrease to the lowest value and then gradually increase. The optimal truncation level J_opt selected by the perturbation process was 23 for N₀ = 150. For N₀ = 180, 120, 90, and 60, the mean BIC_ROI and BIC_sparse series exhibit similar patterns as the degree 150 series (Supplementary Fig. 5). Different numbers J correspond to the lowest BICs for different degrees N₀ are listed in Table 1.

Table 1 Different numbers J correspond to the lowest BICs for different degrees N₀ from the perturbed selection process.

At the completion of the first step, we found the optimal truncation levels J_opt for the different degrees N₀. In the second step, we focused on the selection of the maximum SH expansion degree based on the behavior of the BIC_ROI and BIC_sparse series. As shown in Table 1 and Supplementary Fig. 5, the numbers J corresponding to the lowest value for the BIC_ROI series and the values for the BIC_sparse series are very close to each other or identical for N₀ = 60, 90, 120, and 150. For N₀ = 180, there is a big discrepancy between these two values of J (28 and 17). Therefore, the maximum expansion degree N was selected to be the one with the highest maximum SH degree N₀ for which the numbers J_opt determined by the BIC_ROI and BIC_sparse series were consistent. In our case, the maximum SH expansion degree N was chosen to be 150.

For the maximum degree N = 150 (d150), using the true gravity signals at the repeated survey gravity stations (Supplementary Fig. 6a) and J_opt = 23, the recovered d150 gravity field, ${g}_{{{{{{\rm{model}}}}}}}^{{{{{{\rm{d}}}}}}150}$, is shown in Supplementary Fig. 6b. The residual between ${g}_{{{{{{\rm{model}}}}}}}^{{{{{{\rm{d}}}}}}150}$ and the original d150 field, ${g}_{{{{{{\rm{true}}}}}}}^{{{{{{\rm{d}}}}}}150}$ (Supplementary Fig. 3d), is shown in Supplementary Fig. 6c. The correlation coefficient between ${g}_{{{{{{\rm{model}}}}}}}^{{{{{{\rm{d}}}}}}150}$ and ${g}_{{{{{{\rm{true}}}}}}}^{{{{{{\rm{d}}}}}}150}$ is 0.984. We also computed the variance reduction (VR; see Han et al.⁶⁰ and Rahimi et al.⁶¹): VR = [1 − var(${g}_{{{{{{\rm{true}}}}}}}^{{{{{{\rm{d}}}}}}150}$ − ${g}_{{{{{{\rm{model}}}}}}}^{{{{{{\rm{d}}}}}}150}$)/var(${g}_{{{{{{\rm{true}}}}}}}^{{{{{{\rm{d}}}}}}150}$)] × 100%, where var() is an operator calculating the variance. The derived ${g}_{{{{{{\rm{model}}}}}}}^{{{{{{\rm{d}}}}}}150}$ agrees with ${g}_{{{{{{\rm{true}}}}}}}^{{{{{{\rm{d}}}}}}150}$, yielding a VR of more than 96.7%. These good agreements between ${g}_{{{{{{\rm{model}}}}}}}^{{{{{{\rm{d}}}}}}150}$ and ${g}_{{{{{{\rm{true}}}}}}}^{{{{{{\rm{d}}}}}}150}$ imply that the gravity field recoveries using J_opt = 23 and the d150 true gravity signals are sufficiently accurate in the region of interest.

Broadband leakage mitigation

At this point, we have finished recovering the degree N₀ gravity fields using the synthetic degree N₀ gravity signals at the repeated stations. However, the gravity signals at the Earth’s surface, i.e., the terrestrial gravity observations, contain the spectrum up to an infinite degree (N₀ = $\infty$). The difference between the given recovery degree (N₀) and that contained in the input signals leads to an issue called broadband leakage, which can be described as follows: when recovering a degree N₀ field, the high degree (>N₀) components in the input signal will affect the estimation of the low degree (≤N₀) coefficients. This broadband leakage issue has been reported in previous studies^34,55,57,62. However, few studies have addressed the mitigation of broadband leakage⁵⁷.

We implemented an iterative strategy based on the work of Han et al.⁵⁷ to mitigate broadband leakage. First, the gravity signals of degree 360 (d360) at the repeated survey gravity stations were generated to represent the real terrestrial gravity measurements (Supplementary Fig. 6d) based on the d360 gravity changes from WGHM (Supplementary Fig. 3b). Note that degree 360 is the maximum degree for gravity changes that we can obtain from the WGHM outputs (spatial resolution of 0.5°). Next, we recovered the d150 gravity field using these d360 measurements. The recovered d150 field is shown in Supplementary Fig. 6e, while the residual between the recovered and original d150 fields is shown in Supplementary Fig. 6f. The recovered d150 field using the d360 measurements (Supplementary Fig. 6e) differs notably from the result recovered using the d150 measurements (Supplementary Fig. 6b) in both spatial pattern and strength, especially in the regions around Beijing, Taiyuan-Shijiazhuang-Hengshui, and Zhengzhou. Large discrepancies are also evident in Supplementary Figs. 6c and 6f. The large discrepancies between Supplementary Figs. 6e and 6b, as well as 6c and 6f highlight the existence and effects of broadband leakage.

Next, we removed the high-degree components from the d360 measurements using the Gaussian spatial averaging approach. Assuming that ${\bar{g}}_{j}$ are the Gaussian-weighted average values, we have⁶³

$${\bar{g}}_{j} = \frac{{\sum }_{{{{{{\rm{i}}}}}}=1}^{k}{g}_{{{{{{\rm{i}}}}}}}\cdot {{{{{\rm{W}}}}}}\left({\psi }_{i}\right)}{{\sum }_{{{{{{\rm{i}}}}}} =1}^{k}{{{{{\rm{W}}}}}}\left({\psi }_{i}\right)}\\ W\left(\psi \right) = \frac{2b{e}^{-b\left(1-\cos \psi \right)}}{1-{{{{{{\rm{e}}}}}}}^{-2b}}\\ b = \frac{{{{{{\mathrm{ln}}}}}}\,2}{\left(1-\cos \left(r/R\right)\right)}$$

(3)

where r is the averaging radius, ${\psi }_{i}$ is the spherical distance between the points ${g}_{i}$ and ${g}_{j}$ (0 ≤ ${\psi }_{i}$ ≤ π), $W(\psi )$ are the Gaussian averaging functions, ${g}_{i}$ $(i=1,...,k)$ are the gravity measurements falling in the averaging radius, R is the Earth’s mean radius (6371 km), and b is a dimensionless parameter.

To mitigate the effects caused by the broadband leakage, removing the high degree (>N₀) components from the d360 measurements, as much as possible, is necessary. An iterative process achieves this.

(1) For the averaging radius r, an initial value of 80 km is set. Then, the Gaussian-weighted average $\bar{g}$ is derived according to Eq. 3.

(2) The d150 gravity field is recovered based on $\bar{g}$, and the result is denoted ${g}_{{{{{{\rm{model}}}}}}}^{({{{{{\rm{i}}}}}})}$, where i indicates the iteration number.

(3) The root mean square (RMS) difference between ${g}_{{{{{{\rm{model}}}}}}}^{({{{{{\rm{i}}}}}})}$ and the original field ${g}_{{{{{{\rm{true}}}}}}}^{{{{{{\rm{d}}}}}}150}$ is calculated.

(4) The next iteration begins by increasing the averaging radius r by 1 km.

The RMS difference reaches the smallest value when r equals 99 km (after 20 iterations). After the 99-km Gaussian spatial averaging is applied to the d360 measurements, the derived $\bar{g}$ values (Supplementary Fig. 6g) are consistent with the d150 measurements (Supplementary Fig. 6a). We then recovered the d150 gravity field again. The final d150 field is shown in Supplementary Fig. 6h, while the residual between the recovered and original d150 fields is shown in Supplementary Fig. 6i. The recovered d150 field using the d360 measurements after the Gaussian spatial average agrees well with the result recovered using the d150 measurements (Supplementary Fig. 6b), suggesting that the broadband leakage has effectively been mitigated. Despite the limitation of using the d360 synthetic gravity changes to represent the real time variable gravity measurement, the strategy proposed here represents a major advance in mitigating broadband leakage in Slepian analyses.

Comparison of degree 150 gravity time series and independent absolute gravity measurements

To assess the derived d150 gravity field solutions, we compared the d150 gravity changes with in situ absolute gravity (AG) measurements at six gravity stations in the region of interest. The locations of these stations are shown in Fig. 5a. The AG measurements at station 1 were collected by our group using FG5X-249 and A10-034 gravimeters. The AG measurements were made with FG5/5X and A10 gravimeters for the rest of the stations and were provided by the National Earthquake Data Center (http://data.earthquake.cn). In general, no less than 25 h of repeated measurements were collected at each station, with 25 and 50 accepted sets for the FG5/5X and A10 gravimeters in each measurement. Each set includes 100 free-fall drops. The final gravity values were transferred from the instrumental height to the ground. The accuracy of these AG measurements ranges from 0.60 μGal to 4.09 μGal. None of the AG measurements were included in our d150 gravity field recovery procedure or the previous network adjustment procedure.

In general, the d150 gravity changes at the six stations corresponded well with the in situ AG measurements (Fig. 7). Large discrepancies (approximately 10 μGal) were found in 2010 at stations 4 and 5 (Fig. 7d, e). Because of the limited and sparse in situ AG measurements, the assessment shown here may be subject to large uncertainties. Moreover, the spatial resolutions and spectral coverage differ between the AG measurements and the d150 gravity changes. (1) The AG measurements at a single station primarily represent local mass changes within a limited area and depend on the local hydrogeological characteristics. In contrast, the d150 gravity changes reflect mass changes within an area determined by the spatial resolution of the degree 150 fields (180°/150 = 1.2° in half-wavelength). (2) The AG measurements contain the spectrum up to an infinite degree, while the d150 gravity changes are truncated at degree 150. These spatial resolutions and spectral coverage differences may result in the mismatch between the AG measurements and the d150 time series in 2010 at stations 4 and 5. The assessment shows that the gravity changes from our d150 fields agree reasonably well with the in situ AG measurements. A thorough assessment of IGP-NorthChina2022TG, however, requires terrestrial gravity measurements at additional stations with longer temporal coverages.

Fig. 7: Gravity variations at the gravity stations derived from degree 150 IGP-NorthChina2022TG and in situ absolute gravity (AG) measurements.

Comparisons at a station 1; b station 2; c station 3; d station 4; e station 5; and f station 6. The locations of these stations are denoted by yellow crosses in Fig. 5a. Our gravity field recovery procedure does not include these in situ AG measurements. The red error bars indicate the accuracy of absolute gravity measurements. Offsets are added to the AG-derived results compared with the degree 150 estimates.

Gravity changes from GRACE

During the period of the terrestrial gravity measurements, the September 2016 GRACE data were missing. We used the August 2016 GRACE data instead. As a common postprocessing practice, we replaced the degree-two order-zero coefficients with the more reliable estimates from satellite laser ranging⁶⁴ and added back the degree-one coefficients⁶⁵. After the above corrections, the gravity changes were computed using Slepian functions from the SH coefficients^66,67, with the truncation number J_opt = 3 (same as the d60 terrestrial case). When estimating mass changes using GRACE SH coefficients, a filtering procedure is often necessary to suppress the errors in these coefficients. However, the Slepian method can extract information from GRACE data without further filtering or smoothing^66,67. Next, we estimated the rates of the gravity changes derived using the Slepian functions. We adopted an ensemble mean of the rates from the three data centers to reduce the uncertainties in the GRACE data⁶⁸. Because the ensemble GRACE rates were derived using Slepian functions (Fig. 2d), we also obtained ensemble GRACE rates using two popular DDK⁶⁹ and Gaussian smoothing⁷⁰ filters. The ensemble rates derived using DDK5 and 200-km Gaussian smoothing are shown in Figs. 2e and 2f, respectively. Other DDK and Gaussian smoothing alternatives yielded similar spatial patterns in the region of interest (see Supplementary Note 3, Supplementary Figs. 10 and 11).

Gravity effects induced by vertical displacements

Based on the GPS-derived vertical displacements in North China, we accounted for the gravity effects induced by vertical displacements using the following relation:^39,71

$$\frac{{{{{{{\rm{d}}}}}}g}_{{{{{{\rm{res}}}}}}}^{{{{{{\rm{degree}}}}}}}}{{{{{{\rm{d}}}}}}t}=\frac{{{{{{\rm{d}}}}}}{g}_{{{{{{\rm{TG}}}}}}}^{{{{{{\rm{degree}}}}}}}}{{{{{{\rm{d}}}}}}t}-{\gamma }_{g}\frac{{{{{{{\rm{d}}}}}}h}_{{{{{{\rm{GPS}}}}}}}^{{{{{{\rm{degree}}}}}}}}{{{{{{\rm{d}}}}}}t}$$

(4)

where TG stands for terrestrial gravity, degree indicates the maximum expansion degree (d60 or d150), res indicates the residual term, and ${\gamma }_{g}$ is the vertical gravity gradient. We compared three gravity gradients in North China, namely the free-air gradient (${\gamma }_{g}$ = –3.086 μGal cm^–1), the Bouguer gradient (${\gamma }_{g}$ = –1.967 μGal cm^–1), and the Bouguer-corrected free-air gradient of a spherical body (BCFAG_spherical). Compared to the free-air and Bouguer gradients, the BCFAG_spherical is more realistic in volcanic context and thus widely adopted to address gravity changes of a magma reservoir modeled as a spherical body^33,40,41. Taking into consideration the mass accompanying the deformation, the BCFAG_spherical depicts the changes in gravity that result from an elevation change^40,41

$${\gamma }_{g}={{{{{\rm{free}}}}}}{\mbox{-}}{{{{{\rm{air\; gradient}}}}}}+\frac{4{{{{{\rm{\pi }}}}}}G{\rho }_{0}}{3}\times {10}^{6}\,{{{{{\rm{\mu }}}}}}{{{{{\rm{Gal}}}}}}\,{{{{{{\rm{cm}}}}}}}^{{-}1}$$

(5)

where G is the gravitational constant (6.674 × 10^-11 N m² kg^–2) and ${\rho }_{0}$ is the density of the source in units of kg m^-3.

Since the land subsidence in the NCP has a strong relationship with groundwater overexploitation⁵, we assume a spherical body of groundwater by analogy with the volcanic case when applying the BCFAG_spherical correction. The BCFAG_spherical results in ${\gamma }_{g}$ = –2.806 μGal cm^–1 with groundwater density ${\rho }_{0}$ = 1000 kg m^-3.