Improving the Resolution of GRACE/InSAR Groundwater Storage Estimations Using a New Subsidence Feature Weighted Combination Scheme

: GRACE observations and land subsidence data derived from InSAR both assess groundwater storage changes. However, GRACE data at local scales are restricted by the coarser spatial resolution of satellite systems, and inversion of Groundwater Storage Anomalies ( GWSA ) by InSAR requires extensive and unavailable lithological data. Here, we propose a New Subsidence Feature Weighted Combination (NSFWC) scheme to enhance the spatial resolution of GRACE-derived GWSA from 0.5 ◦ to 0.05 ◦ . This method can not only retain the spatial distribution of groundwater changes but also reﬂect local details related to surface subsidence. A case study was executed to evaluate the performance of the NSFWC scheme in the Beijing Plain, which has seriously overexploited groundwater. Results showed that the simulated GWSA were consistent with in situ measurements in most regions, with a correlation coefﬁcient of 0.85 and an RMSE of 4.41 mm/year. Additionally, there were 22 overexploited wells in the Beijing Plain, although groundwater levels generally recovered after the South to North Water Diversion Project. Simultaneously, four cones of depression were detected by the InSAR technology, where the maximum cumulative subsidence and subsidence rate achieved − 198.52 mm and − 53.09 mm/year, respectively. This paper provides data support and technical guarantees for small-scale groundwater resources management.


Introduction
Many semi-arid regions suffer from intensive groundwater depletion that exceeds natural replenishment. Excessive use of groundwater resources increases the risk of decreasing aquifers' pore water levels. The changes of effective stress result in surface deformation because of the compaction of cohesive soil. The consequences of subsidence may reduce flood discharge efficiency and damage urban infrastructures such as buildings, highways, and bridges [1]. Additionally, the inelastic sediments caused by aquifer-system compaction lead to the inability of aquifer storage and unsustainable water management [2]. The problem of subsidence induced by groundwater overexploitation has become increasingly prominent in the North China Plain, such as in Zhengzhou [3], Beijing [4], Tianjin [5], Dezhou [6], and Changzhou [7]. An adequate measuring tool for water authorities is necessary to solve supply and demand problems.
Currently, traditional groundwater monitoring approaches comprise in situ wells and groundwater budgets. However, field points are sparse, and monitoring is timeconsuming and labor-intensive. Thus, drilling observations are insufficient for capturing the ilating measurements into existing process-based models. Statistical downscaling methods attempt to derive empirical relationships between GRACE observations and smaller-scale quantities of interest. Previous works, such as Zhong et al. [29], Arshad et al. [30,31], Ali et al. [32,33], and Noor et al. [34], have reported the results of downscaling in the hydrological field. For example, Arshad et al. (2021) made an attempt to develop an integrated downscaling and calibration framework to generate high-resolution gridded precipitation data. The result reported that a Mixed Geographically Weighted Regression (MGWR) model, capable of dealing with fixed and spatially varied environmental variables, was proposed to downscale the original TRMM precipitation data from 0.25 • to 1 km [30]. Ali et al. (2023) employed machine learning models, such as extreme gradient boosting and artificial neural networks, to downscale GRACE-TWSA from 1 • to 0.25 • [33]. Several attempts have been made to combine both approaches to strengthen the capabilities of In-SAR and GRACE [35]. For instance, Castellazzi et al. (2016) expounded on the benefits and challenges of integrating GRACE and InSAR and reported the theoretical and conceptual basis of the two approaches [13]. Subsequently, Castellazzi et al. (2018) found that ground displacement could be regarded as a prior spatial map of groundwater storage loss [36]. Shang et al. (2019) designed a downscaling model using the relevance vector machine, inversing the overexploitation-induced ground subsidence [37]. Liu et al. (2019) observed the displacement resulting from groundwater pumping in the southern Central Valley using Sentinel-1 data. Results showed a good correlation between GRACE-derived GWSA and sedimentation records [34]. Massoud et al. (2021) measured the groundwater changes and estimated the land movement based on the high-resolution InSAR algorithm [38]. Bai et al. (2022) investigated the GWSA in Cangzhou in the past decades by coupling multiple SAR images and hydraulic head measurements [7]. However, an effective method has not yet been proposed to improve the spatial resolution by combining both data in the Beijing Plain.
According to previous research, land subsidence rates were extremely interrelated to the declining trend of groundwater levels over the overexploitation centers [13,39]. However, as one of the effective methods to detect groundwater changes, the coarse resolution of GRACE restricted its application at a local scale. Moreover, the spatial distribution of groundwater changes monitored by GRACE was usually inconsistent with official data due to leakage errors. Little research has been published on the GWSA of aquifer compaction by combining InSAR-derived deformation in spatial patterns. Nevertheless, the spatial downscaling performance of machine learning will be affected by the input variables, applied algorithms, and data period. This paper proposed the New Subsidence Feature Weighted Combination (NSFWC) scheme to investigate GRACE signals at the resolution of 0.05 • . This method reflects the physical processes with high spatial resolution, and it is suitable to be popularized because it does not need a large amount of auxiliary data. The Beijing Plain, where the groundwater is seriously exhausted, was utilized to confirm the effectiveness of the NSFWC scheme. Section 1 introduces the research background and significance. Section 2 describes the hydrogeological conditions of the Beijing Plain, the data processing methods of in situ observations, InSAR, and GRACE, as well as the theory of the feature-weighted combination scheme. Section 3 illustrates the results of in situ observations, land subsidence, groundwater storage anomaly, and simulated GWSA from GRACE/InSAR. The impact factors of subsidence and simulated uncertainties of vertical displacement are discussed in Section 4. The conclusions are revealed in Section 5.
The main goal of this study is to downscale the spatial resolution of GRACE-derived GWSA from 0.5 • to 0.05 • by merging InSAR data. The sub-objectives include (1) obtaining field data to display the true groundwater level changes, (2) studying land subsidence using time-series InSAR technology, (3) retrieving the GWSA using GRACE and GLDAS data, and (4) designing an optimum fusion model based on (2) and (3), as well as comparison with field data.

Study Area
Beijing belongs to the Haihe River Basin along with Tianjin and Hebei. It is located at 115.5 E~117.5 E and 39 N~41 N, as shown in Figure 1. The Beijing Plain consists of 13 districts, including the 8 districts in the old city, as well as portions of Pinggu, Tongzhou, Miyun, Huairou, and Shunyi. The spatiotemporal distribution of precipitation is uneven, with relatively rainy mountains and rainless plains. It suffered drought for five consecutive years beginning in 1999, which resulted in a main contradiction between water supply and water use. Afterward, the aquifer system was perennially overexploited to ensure water supply. The average groundwater depth in the Beijing Plain declined from 11.9 m in 1998 to 24.3 m in 2012 [40]. Rapid ground deformation was caused by groundwater pumping in past decades.
using time-series InSAR technology, (3) retrieving the GWSA using GRACE and GLDA data, and (4) designing an optimum fusion model based on (2) and (3), as well as compar ison with field data.

Study Area
Beijing belongs to the Haihe River Basin along with Tianjin and Hebei. It is located at 115.5 E~117.5 E and 39 N~41 N, as shown in Figure 1. The Beijing Plain consists of 1 districts, including the 8 districts in the old city, as well as portions of Pinggu, Tongzhou Miyun, Huairou, and Shunyi. The spatiotemporal distribution of precipitation is uneven with relatively rainy mountains and rainless plains. It suffered drought for five consecu tive years beginning in 1999, which resulted in a main contradiction between water suppl and water use. Afterward, the aquifer system was perennially overexploited to ensur water supply. The average groundwater depth in the Beijing Plain declined from 11.9 m in 1998 to 24.3 m in 2012 [40]. Rapid ground deformation was caused by groundwate pumping in past decades. The compressible materials of aquifers form an overlying stratum, and the compac tion thickness increases from west to east. According to the Quaternary sedimentary pat tern, aquifer structure, and groundwater consumption, aquifer systems in Beijing are di vided into four layer groups [41]. The lithology of the phreatic aquifer mainly consists o silt, while that of confined aquifers is composed of clay, silt, and sand. Remarkably, th third confined aquifer is the major extraction layer for living. The groundwater level o the aquitard generally appears as elastic-plastic changes, which is the primary contrib uting layer of long-term subsidence [37]. The operation of the SNWD project is expected to restore the groundwater level through ecological replenishment and irrigated infiltra tion. The government has strictly limited groundwater pumping in recent years; however some regions are still overexploited. Therefore, detailed groundwater monitoring on a lo cal scale is necessary for comprehending groundwater resources. The compressible materials of aquifers form an overlying stratum, and the compaction thickness increases from west to east. According to the Quaternary sedimentary pattern, aquifer structure, and groundwater consumption, aquifer systems in Beijing are divided into four layer groups [41]. The lithology of the phreatic aquifer mainly consists of silt, while that of confined aquifers is composed of clay, silt, and sand. Remarkably, the third confined aquifer is the major extraction layer for living. The groundwater level of the aquitard generally appears as elastic-plastic changes, which is the primary contributing layer of long-term subsidence [37]. The operation of the SNWD project is expected to restore the groundwater level through ecological replenishment and irrigated infiltration. The government has strictly limited groundwater pumping in recent years; however, some regions are still overexploited. Therefore, detailed groundwater monitoring on a local scale is necessary for comprehending groundwater resources.

SAR Data
In total, 39 Sentinel-1A single-look complex images were obtained from the European Space Agency [42]. The SAR datasets were collected at intervals of approximately one month, and the missing values for March and September 2019 were filled using linear interpolation. The Sentinel-1A satellite includes 4 scanning modes, including Stripmap (SM), Interferometric Wide swath (IW), Extra-Wide swath (EW), and Wave (WV). The incident angle of Sentinel-1A data is approximately 39.5 • in the interference width swath mode. The terrain phase was removed with the 30 m SRTM DEM. At the same time, POD precise orbit data were downloaded as the auxiliary data for SAR pre-processing. We summarized the parameters of Sentinel-1A, as shown in Table 1.

GRACE and GLDAS Data
Two new versions of GRACE regularized mascon solutions from May 2017 to July 2020 were considered: the Jet Propulsion Laboratory RL06 v02 (JPL) [43] and the Center for Space Research (CSR) [44]. Other components needed to be removed from Terrestrial Water Storage Anomalies (TWSA) when GRACE was applied to retrieve GWSA. Soil Moisture Storage Anomalies (SMSA) and Surface Water Storage Anomalies (SWSA) were both calculated using the Global Land Data Assimilation System (GLDAS) NOAH land surface model. For the sake of spatial compatibility of different datasets, the resolution was unified to 0.05 • × 0.05 • . All components derived using GRACE and NOAH are required to be interpolated and anomaly processed.

Validation Data
The groundwater level data originate from the Statistical Yearbook of Groundwater Level in China's Geological Environment. There are 102 monitoring wells distributed in the Beijing Plain. Local governments manage the volume of groundwater exploitation according to aquifer availabilities. The groundwater level in water-receiving regions gradually recovered after the operation of SNWD in December 2014. The pumping amount of different aquifers should be reasonably arranged according to the residuals of recharge and discharge fluxes. GPS stations were collected from the Crustal Movement Observation Network of China. Two GPS sites (BJSH and BJFS) in Beijing were considered to confirm the reliability of InSAR. However, antenna replacement or earthquake may induce a jump mutation of sequences at a certain time, resulting in the uneven sampling of GPS stations. It is essential to carry out pre-processing before analyzing original GPS data, such as anomaly detection, temporal interpolation, and leap correction [45].

Method
The purpose of this study was to enhance the resolution of GRACE-derived GWSA from 0.5 • to 0.05 • using the NSFWC scheme. Since InSAR cannot directly provide quantitative groundwater changes, it is difficult to use it alone. Firstly, the distribution of land deformations was detected by InSAR technology in pumping regions. Secondly, groundwater storage anomalies were derived from GRACE data. Thirdly, as the prior spatial knowledge, deformations were combined with groundwater storage anomalies into a new set of estimations. Indeed, GRACE-based GWSA with a higher resolution was acquired, which could satisfy the need for water resources management. Considering the quality and availability of data, as well as the geological conditions of the Beijing Plain, the time span during 2017-2020 was selected as the study period. Figure 2 shows the data processing flowchart.
The purpose of this study was to enhance the resolution of GRACE-derived GWSA from 0.5° to 0.05° using the NSFWC scheme. Since InSAR cannot directly provide quantitative groundwater changes, it is difficult to use it alone. Firstly, the distribution of land deformations was detected by InSAR technology in pumping regions. Secondly, groundwater storage anomalies were derived from GRACE data. Thirdly, as the prior spatial knowledge, deformations were combined with groundwater storage anomalies into a new set of estimations. Indeed, GRACE-based GWSA with a higher resolution was acquired, which could satisfy the need for water resources management. Considering the quality and availability of data, as well as the geological conditions of the Beijing Plain, the time span during 2017-2020 was selected as the study period. Figure 2 shows the data processing flowchart.

InSAR Measurements
The improved multi-temporal InSAR approach based on SBAS-InSAR and PS-InSAR was utilized to survey the spatiotemporal assessment of land deformation, and the model was constructed as follows [46].

InSAR Measurements
The improved multi-temporal InSAR approach based on SBAS-InSAR and PS-InSAR was utilized to survey the spatiotemporal assessment of land deformation, and the model was constructed as follows [46]. where φ int refers to the interference phase; φ re f stands for the reference ellipsoid phase; φ dem suggests the topographical phase, which refers to the relative ground motion during the two shooting dates. φ atm is the atmospheric delay phase, and φ res means the residual phase. B 0 and B 0 ⊥ represent the parallel and vertical projection components of the space baseline in the line of sight (LOS) direction, respectively. θ means incidence angle, λ indicates the radar wavelength, R is the radar slant distance, ∆R signifies the displacement in LOS, and h is ground elevation relative to the reference ellipsoid. Note that atmospheric delay signals are not eliminated by the high-pass filter because it contains seasonal information. Therefore, the time-series deformation may contain contributions of the atmosphere, which is low due to the flat topography in the Beijing Plain.
Long-term displacement based on the InSAR method consists of the following crucial steps. Firstly, SAR images were processed using differential interferometry [47]. The SAR image on May 2017 was selected as the super master (SM) image, and the others images were re-sampled and registered to SM. During the interference process, the multilook processing of 10 × 2 was used to improve the signal/noise ratio. Secondly, DEM data were employed to reduce the φ dem from the interference phase after interferograms and amplitude values were generated. Targets with high coherence and stable scattering properties were taken as persistent scatterers, such as road edges, bridge bodies, exposed rocks, and house roofs. The coherence of interferograms was evaluated by the coherence coefficient threshold (0 < γ < 1) combining intensity and coherent information, which was set as 0.3 in this paper. Thirdly, a three-dimensional phase unwrapping method was utilized to regain the true values. Fourthly, the Goldstein filtering algorithm was used to suppress noise, and then the linear phase slope of each interferogram was subtracted to eliminate orbit errors and atmospheric impacts. Finally, time-series LOS deformations were ultimately retrieved with the least square method and geocoded to the geographic coordinate system. The deformations measured by InSAR were calculated relative to the reference point, which was Yuyuantan Park, almost free of compression. Since the reference point had stable backscattering characteristics, its deformation was assumed to be zero. Owing to the large coverage area of the Beijing Plain, a total of 11 bursts were processed with the aid of the reference point to unify the benchmark. Generally speaking, ground deformation in Beijing is usually caused by groundwater extraction, so horizontal displacement caused by crustal movement is very small compared with vertical displacement. In addition, the radar wave is the most sensitive to vertical displacement, followed by east-west displacement, and the least sensitive to north-south motion. Therefore, horizontal deformations (the east-west and north-south directions) were negligible. According to the incident angle and radar parameters, LOS deformations were converted into perpendicular subsidence by the trigonometric function. The calculation method was as follows [48,49].
where D LOS is the LOS deformation; D v is the perpendicular movement; D n and D e indicate the horizontal movements; θ, α, and ϕ represent radar incidence angle, heading angle, and azimuth angle, respectively; δ LOS is radar noise.

GRACE Solutions
Generally, TWSA contains five components, including GWSA, SMSA, SWSA, Canopy Water Storage Anomalies (CWSA), and Snow Water Storage Anomalies (SnWSA). For the Beijing Plain, the SnWSA are relatively small, so they can be ignored [50]. Research [51] implies that changes in water resources caused by plants are about 5 mm, far less than GRACE uncertainties (2 cm). Hence, canopy water storage is negligible, too [50]. Under these hypotheses, the distribution maps of GWSA were obtained by subtracting Soil Mois-ture Storage Anomalies (SMSA) and Surface Water Storage Anomalies (SWSA) from TWSA under the same resolution, which is expressed by the following equation [52,53].
where TWSA k is the terrestrial water storage anomalies of the kth grid from GRACE; SMSA k is the NOAH-derived soil moisture storage anomalies of the kth grid; SWSA k is the NOAH-derived surface water storage anomalies of the kth grid. The Mann-Kendall (MK) [54] test is utilized to assess the significance level of hydrological data. It indicates significant trends with 95% confidence when the absolute value of Z overtakes 1.96. On the contrary, non-significant trends are indicated with grey dots in the spatial distribution map (p ≥ 0.05). Additionally, singular spectrum analysis (SSA) was applied to decompose original data into three components, including long-term trend, periodic season, and residual error [55].

New Subsidence Feature Weighted Combination Scheme
Deformations derived from InSAR provided high-density and evenly distributed pixels, including uplift zones and subsidence zones. The inversion process of the NSFWC scheme depended on finding a group of pixels with a declining deformation rate. The deformation map was taken as a priori knowledge to process GWSA at the grid scale. After the sedimentation pixels were extracted, InSAR-based subsidence and GRACE-based GWSA were converted into new values at the unified resolution to reduce the simulation errors. The NSFWC method was verified by the absolute errors, and the processing workflows were mainly divided into the following five steps.
The pending weight coefficient was obtained by the corresponding minimum kriging variance on the premise of unbiased, which was determined by the following equation.
where σ i stands for the pending weight coefficient; Z * (λ 0 , ϕ 0 , h 0 ) means deformation results based on the kriging interpolation method; Second, 16,353 pixels were generated in Beijing Plain, which was re-sampled into 406 grids. Then, sinking grids were extracted as the data source for calculating weight in the next step [57].
where D sub is the rate value of decreasing grid; D v is the trend of vertical displacement; N is the total grid number, N = 406. Third, the sample data can be processed into dimensionless data in the range of [0, 1] according to the Max-Min normalization algorithm [58]: However, the Max-Min normalization algorithm is easily affected by extreme values. We took the average subsidence rate at the four funnels as the minimum value, and then an improved normalization algorithm was obtained. The subsidence rate was inversely proportional to the weight in the new algorithm. The improved normalized data was assigned to GRACE as the weight of the NSFWC scheme.
where W k is the normalization result of subsidence grids; D A v , D B v , D C v , and D D v are the vertical displacements at four cones of depression (A, B, C, and D) respectively; D k sub denotes the subsidence rate of the kth grid; D max sub means the maximum subsidence rate; and n is the total grid number of depression funnels.
Forth, GRACE-derived GWSA were re-estimated based on the NSFWC scheme. The new algorithm added the characteristics of surface deformation caused by groundwater exploitation.
where GWSA k M represents the GWSA estimated by the feature combined scheme; GWSA k denotes GRACE-based groundwater storage anomaly of the kth grid; k denotes the number of subsidence grids in the study area, k = 1, 2, . . . , 258.
Fifth, The GWSA estimated from groundwater level changes ( GWL) were regarded as reference values, which were applied to affirm the reliability of combined GWSA based on the NSFWC scheme [59].
where GWSA k O suggests the GWSA estimated by groundwater level changes of the kth grid, ∆GW L k is the groundwater level changes of the kth grid, and Sy is the average specific yield in the phreatic aquifer [60].
Ultimately, the simulation results were evaluated by absolute error (e k ) [57], Pearson correlation coefficient (r) [61], and root mean square error (RMSE) [62]: where e k denotes the absolute error at the kth grid, GWSA k M indicates the simulation result at the kth grid, GWSA k O means the in situ data at the kth grid, GWSA M is the average of simulated results, GWSA O is the average of in situ data. m is the number of settlement grids. 2015, the operation of SNWD has relieved water scarcity in water-receiving regions. The water levels in most wells have recovered to a certain degree, and the transformation of the ∆GW L trend also changes from negative to positive (Figure 4). The average trends of groundwater depth in three periods are calculated by the linear fitting method at the rate of −0.91 m/year, 0.16 m/year, and 0.03 m/year, respectively. However, there are still some monitoring wells whose trends continue to decline at the moment, and a total of 22 such points have been detected. Ten points at downward trends are selected to display the temporal GWL changes in the third phase, as shown in Figure 5, which are marked with red triangles in Figure 1. shows decreased trends due to continuous groundwater over exploitation. Then, it becomes gradually stable owing to the increase in rainfall durin 2012-2014. Since 2015, the operation of SNWD has relieved water scarcity in water-receiv ing regions. The water levels in most wells have recovered to a certain degree, and th transformation of the GWL Δ trend also changes from negative to positive (Figure 4). Th average trends of groundwater depth in three periods are calculated by the linear fittin method at the rate of −0.91 m/year, 0.16 m/year, and 0.03 m/year, respectively. However there are still some monitoring wells whose trends continue to decline at the moment, and a total of 22 such points have been detected. Ten points at downward trends are selected to display the temporal GWL changes in the third phase, as shown in Figure 5, which ar marked with red triangles in Figure 1.     Figure 6a reflects the cumulative deformation of the Beijing Plain from 2017 to 2020. For subsequent combined inversion, the vector map is converted into a grid map at a resolution of 0.05 • , as shown in Figure 6b. It is found that major regions are stable without evident deformation. In terms of spatial distribution, deformations are mainly located in the border zone among Chaoyang, Tongzhou, Changping, and Haidian, forming four depression funnels. In terms of time-series changes, the first image (20 May 2017) is taken as reference values to calculate the cumulative displacement. It reaches up to −183.70~−92.63 mm at the four funnels with a maximum subsidence rate of −53.09 mm/year (Figure 7b). The findings in this paper are similar to the previous studies on the NCP, which is the most serious subsidence region in China [63,64]. However, many previous works focused on extracting Beijing's deformation using ENVISAT data from 2007 to 2010, during which the serious subsidence speed reached 150 mm/year [65,66]. Owing to the SNWD replenishment and the strict groundwater pumping policy, ground subsidence in most urban areas was alleviated in 2017-2020.    Figure 6a reflects the cumulative deformation of the Beijing Plain from 2017 to 2020. For subsequent combined inversion, the vector map is converted into a grid map at a resolution of 0.05°, as shown in Figure 6b. It is found that major regions are stable without evident deformation. In terms of spatial distribution, deformations are mainly located in the border zone among Chaoyang, Tongzhou, Changping, and Haidian, forming four depression funnels. In terms of time-series changes, the first image (20 May 2017) is taken as reference values to calculate the cumulative displacement. It reaches up to −183.70~−92.63 mm at the four funnels with a maximum subsidence rate of −53.09 mm/year (Figure 7b). The findings in this paper are similar to the previous studies on the NCP, which is the most serious subsidence region in China [63,64]. However, many previous works focused on extracting Beijing's deformation using ENVISAT data from 2007 to 2010, during which the serious subsidence speed reached 150 mm/year [65,66]. Owing to the SNWD replenishment and the strict groundwater pumping policy, ground subsidence in most urban areas was alleviated in 2017-2020.     Figure 8 presents 258 sinking grids estimated by InSAR, and the blank grids ex the uplift zone. Four monitoring wells (P17, P66, P68, and P86) located at the descend are selected to analyze the relationship between △GWL and subsidence. It is foun trends of InSAR subsidence are consistent with △GWL, JPL-GWSA, CSR-GWSA CLSM-GWSA by comparison (Figure 9a-d). Note that GWSA retrieved with JPL and stands for the annual average, and the time span of data availability is from June 2 December 2020. Additionally, land deformation is delayed by 1 ~ 6 months compared △GWL. To reasonably analyze the relationship between land subsidence and groun ter changes, it is necessary to pre-treat △GWL. The benchmark values (20 May 201 subtracted from the original records, and the trend component is decomposed by th method before comparison. Figure 9e-h show the correlation among their trend co nents of InSAR subsidence, CLSM-GWSA, and △GWL at the four monitoring well correlation coefficients are shown in Table 2, which are above 0.85 during 2017-20 other words, subsidence rates rise with the decrease in groundwater level; meanw subsidence rates also slow down when the groundwater level rises.  Figure 8 presents 258 sinking grids estimated by InSAR, and the blank grids express the uplift zone. Four monitoring wells (P17, P66, P68, and P86) located at the descend zone are selected to analyze the relationship between GWL and subsidence. It is found that trends of InSAR subsidence are consistent with GWL, JPL-GWSA, CSR-GWSA, and CLSM-GWSA by comparison (Figure 9a-d). Note that GWSA retrieved with JPL and CSR stands for the annual average, and the time span of data availability is from June 2018 to December 2020. Additionally, land deformation is delayed by 1~6 months compared with GWL. To reasonably analyze the relationship between land subsidence and groundwater changes, it is necessary to pre-treat GWL. The benchmark values (20 May 2017) are subtracted from the original records, and the trend component is decomposed by the STL method before comparison. Figure 9e-h show the correlation among their trend components of InSAR subsidence, CLSM-GWSA, and GWL at the four monitoring wells. The correlation coefficients are shown in Table 2, which are above 0.85 during 2017-2020. In other words, subsidence rates rise with the decrease in groundwater level; meanwhile, subsidence rates also slow down when the groundwater level rises.        Figure 10 shows the spatial distribution of annual trends, including JPL-TWSA, CSR-TWSA, NOAH-SMSA, and NOAH-SWSA. Figure 10a-d display the spatial trend maps of four components under the original resolution, with 0.5 • for JPL and 0.25 • for the other three components. Then, the trend maps of four components are interpolated to 0.05 • using the cubic spline interpolation method, as shown in Figure 10e-h. The spatial distribution of trends after interpolation is almost the same as before; thus, the errors caused by interpolation are not considered. On the one hand, two GRACE solutions make use of different de-striping strategies, resulting in different effects of de-striped residuals on them. That is the reason why TWSA estimated with JPL is inconsistent with that of CSR (Figure 10a,b). On the other hand, GRACE leakage errors may be caused by data solutions such as spherical harmonic coefficients truncation, de-striping strategies, and Gaussian smoothing filtering. At present, it fails to determine the appropriate GRACE solutions for a specific region. Therefore, three available GWSA derived from JPL, CSR, and CLSM are employed for comparative analysis. them. That is the reason why TWSA estimated with JPL is inconsistent with that of CSR (Figure 10a,b). On the other hand, GRACE leakage errors may be caused by data solutions such as spherical harmonic coefficients truncation, de-striping strategies, and Gaussian smoothing filtering. At present, it fails to determine the appropriate GRACE solutions for a specific region. Therefore, three available GWSA derived from JPL, CSR, and CLSM are employed for comparative analysis. According to relevant studies in Cangzhou, Dezhou, and Langfang, the empirical relationship between subsidence and △GWL generally presents exponential, quadratic, or linear functions [67]. Based on the prior map of land subsidence, GWSA in the Beijing Plain is reconstructed using the NSFWC scheme, which improves the resolution of GRACE-derived GWSA. Figure 11a,b show GWSA trends obtained from two GRACE solutions at the resolution of 0.5° during 2018-2020, according to Equation (1). Figure 11c displays CLSM-GWSA trends during 2017-2020, which are matched with the downloaded SAR images. Three GWSA trends are downscaled to 0.05° based on the NSFWC scheme, as seen in Figure 11d-f. The resolution of simulated GWSA is greatly enhanced compared with the original data, which can meet the needs of water resource management. According to relevant studies in Cangzhou, Dezhou, and Langfang, the empirical relationship between subsidence and GWL generally presents exponential, quadratic, or linear functions [67]. Based on the prior map of land subsidence, GWSA in the Beijing Plain is reconstructed using the NSFWC scheme, which improves the resolution of GRACEderived GWSA. Figure 11a The InSAR/GRACE simulated GWSA performs detailed groundwater storage changes, indicating that InSAR can serve a priori of groundwater depletion map. Furthermore, lithology (aquifer property, thickness, and compressibility) is different over the whole Beijing Plain; there still is a positive correlation between groundwater depletion The InSAR/GRACE simulated GWSA performs detailed groundwater storage changes, indicating that InSAR can serve a priori of groundwater depletion map. Furthermore, lithology (aquifer property, thickness, and compressibility) is different over the whole Beijing Plain; there still is a positive correlation between groundwater depletion rate and subsidence rate. The simulated GWSA from different data sources appear in similar results (Figure 11d-f), which proves the potential of combining GRACE and InSAR. Therefore, the NSFWC scheme proposed in this paper can be applied to future GRACE mascon missions to improve spatial resolution. Figure 12 illustrates the comparison of GWSA, subsidence, and precipitation in terms of long-term trends and seasonal changes. The Beijing Plain is divided into three parts (north, middle, and south), and the average values of each sector are calculated to facilitate the comparative analysis. Note that GWSA is retrieved from GRACE and CLSM; subsidence is derived from InSAR and GPS, and precipitation is sourced from ERA 5. Figure 12a compares the average values of the three kinds of data over the whole Beijing Plain. It is found that long-term trends of InSAR subsidence are consistent with GRACE-GWSA, showing a downward trend, whereas CLSM-GWSA agrees with GPS subsidence without a significant trend. Moreover, the average subsidence rate is higher than the CLSM-GWSA trend from 2015 to 2020. The CLSM model does not reflect the positive impact of SNWD on surface water because of the lack of human activities in the model structure. The discrepancy between subsidence sequences retrieved from InSAR and GPS is that there are few GPS sites in Beijing. There is no doubt that the values of a single station cannot accurately represent displacements at a local scale. while the subsidence oscillates within ± 3 mm. Second, subsidence lags behind GWSA, which is explained by the low perpendicular hydraulic conductivity of aquifer systems and the slow consolidation of soil mass in wet weather [68]. Third, subsidence and precipitation are negatively correlated with a correlation coefficient of −0.45, while that is 0.06 with GWSA, indicating that the vibration is most likely caused by climate change. Because GWSA and subsidence sequences are damaged by noise, errors are produced from the simulated GWSA estimation and seasonal component decomposition. Therefore, we estimate the average values of regional pixels to reduce random noise. It can also be seen from relationships that seasonal components are less relevant than longterm trend components. Figure 13 shows the comparison of InSAR subsidence between individual grids near site BJSH and average values in the northern Beijing Plain. Results reflect that the time series of average values are similar to that of the adjacent grids re-   Figure 12c shows the comparison of seasonal signals, and the results reveal the following three aspects. First, three data consisting of GWSA, settlement, and precipitation all perform obvious seasonal fluctuations with different amplitude. The ranges of precipitation and GWSA are −62.22~87.27 mm and −9.69~13.14 mm, respectively, while the subsidence oscillates within ±3 mm. Second, subsidence lags behind GWSA, which is explained by the low perpendicular hydraulic conductivity of aquifer systems and the slow consolidation of soil mass in wet weather [68]. Third, subsidence and precipitation are negatively correlated with a correlation coefficient of −0.45, while that is 0.06 with GWSA, indicating that the vibration is most likely caused by climate change.

Comparison of GWSA, Land Subsidence, and Precipitation
Because GWSA and subsidence sequences are damaged by noise, errors are produced from the simulated GWSA estimation and seasonal component decomposition. Therefore, we estimate the average values of regional pixels to reduce random noise. It can also be seen from relationships that seasonal components are less relevant than long-term trend components. Figure 13 shows the comparison of InSAR subsidence between individual grids near site BJSH and average values in the northern Beijing Plain. Results reflect that the time series of average values are similar to that of the adjacent grids regardless of the different magnitude. Therefore, it is feasible to estimate subsidence at the plain scale. Because GWSA and subsidence sequences are damaged by noise, errors are produced from the simulated GWSA estimation and seasonal component decomposition. Therefore, we estimate the average values of regional pixels to reduce random noise. It can also be seen from relationships that seasonal components are less relevant than longterm trend components. Figure 13 shows the comparison of InSAR subsidence between individual grids near site BJSH and average values in the northern Beijing Plain. Results reflect that the time series of average values are similar to that of the adjacent grids regardless of the different magnitude. Therefore, it is feasible to estimate subsidence at the plain scale.

Impact of Land Subsidence
Groundwater storage variabilities play a vital role in land subsidence, which is a response to human and climate pressure [15]. On the one hand, the long-term trends of deformation are generally related to aquifer compaction. The specific compression of cohesive soil passes that of sand and gravel under the same exploited conditions. The subsidence unit of aquifers dominated by silt and clay is 2~3 times that of gravel, which is the primary contribution of deformation in Beijing [69]. Tiny particles play a crucial role in preventing perpendicular groundwater flow. The heterogeneity of the displacement rate is induced by the control of early fault events on the sedimentary environment on both sides. This phenomenon triggers different thicknesses of sediment and then leads to the difference in aquifer compression during groundwater exploitation. Occasionally, a large deformation will emerge even if the Quaternary sediment thickness is similar along the same side of the fault. This can be explained by the separation of two different compressible aquifer materials in fault activity. The response of geological media to effective stress may exceed the pre-consolidation stress in several monitoring wells that penetrate the confined aquifer. The time-series subsidence trends result from fine particle rearrangement, which results in the inability of water storage. Concerning a more complete description of ground deformation caused by the compressibility of aquifer systems, the article of Galloway (2007) is beneficial for understanding the mechanism [70].
On the other hand, the elastic change of subsidence is related to the seasonal change of precipitation. Generally, InSAR results are interfered with by various noises, such as atmospheric phase delay [71], which absorbs red noise of higher frequency and produces fluctuations similar to seasonal signals [72]. Thus, the observed subsidence sequence has measurable seasonal fluctuations, and the seasonal signals can be decomposed and compared with precipitation. The singular spectrum analysis method was applied to obtain the seasonal component by eliminating the trend component. Note that the extracted season signals of precipitation may contain noise, such as the outliers caused by drought and excessive rainfall. Generally, there is a delay of less than one year from groundwater level changes to subsidence.
Google Earth shows that there are many vegetable greenhouses, industrial plants, and residential buildings at the cones of depression. As shown in Figure 14a-d, regions with serious subsidence at funnels A, B, C, and D are selected to display their geomorphic characteristics. This indicates that human activity is one of the reasons for ground sedimentation, such as industrial and agricultural production. Overall, these factors consist of groundwater exploitation, human activities, precipitation, and hydrogeology, which jointly affect land subsidence.
tion of ground deformation caused by the compressibility of aquifer systems, the a of Galloway (2007) is beneficial for understanding the mechanism [70].
On the other hand, the elastic change of subsidence is related to the seasonal ch of precipitation. Generally, InSAR results are interfered with by various noises, su atmospheric phase delay [71], which absorbs red noise of higher frequency and prod fluctuations similar to seasonal signals [72]. Thus, the observed subsidence sequenc measurable seasonal fluctuations, and the seasonal signals can be decomposed and pared with precipitation. The singular spectrum analysis method was applied to o the seasonal component by eliminating the trend component. Note that the extracted son signals of precipitation may contain noise, such as the outliers caused by drough excessive rainfall. Generally, there is a delay of less than one year from groundwater changes to subsidence.
Google Earth shows that there are many vegetable greenhouses, industrial pl and residential buildings at the cones of depression. As shown in Figure 14a-d, reg with serious subsidence at funnels A, B, C, and D are selected to display their geomor characteristics. This indicates that human activity is one of the reasons for ground mentation, such as industrial and agricultural production. Overall, these factors cons groundwater exploitation, human activities, precipitation, and hydrogeology, w jointly affect land subsidence.

Uncertainty of Simulated GWSA
The spatial residual map of simulated results and official data is evaluated according to Equation (11). The simulated CLSM-GWSA is utilized for residual analysis to unify the period, as shown in Figure 15. According to the distribution pattern, it can be observed that the majority of errors are below 10 mm/year, and the largest errors happen in the northeast corner and the west of the Beijing Plain. The GWSA signals in these regions fail to explain the mass map derived from InSAR. In addition, subsidence located in the northeast corner and west area cannot be completely attributed to groundwater change, and it may be greatly affected by the hydrological process. It is worth mentioning that there are several separate studies on GWSA [23] and subsidence [37], respectively, over the Beijing Plain, and their reported results are similar to this article. The spatial correlation analysis of the remaining areas was carried out after removing these two abnormal regions. We found that the simulated results were consistent with in situ wells with a correlation coefficient of 0.85 and an RMSE of 4.41 mm/year. The correlations of original GWSA from JPL, CSR, and CLSM data are 0.92, 0.88, and 0.86, respectively. several separate studies on GWSA [23] and subsidence [37], respectively, over th Plain, and their reported results are similar to this article. The spatial correlation of the remaining areas was carried out after removing these two abnormal regi found that the simulated results were consistent with in situ wells with a correla efficient of 0.85 and an RMSE of 4.41 mm/year. The correlations of original GW JPL, CSR, and CLSM data are 0.92, 0.88, and 0.86, respectively. The diversity between simulated GWSA and in situ monitoring wells ca plained from four aspects. At first, GRACE-GWSA does not fully reflect the groundwater depletion due to leakage errors. Second, 36% of grids are not calcu the GWSA trends map at 0.25° original resolution ((406-258)/406 ≈ 0.36), which me several grids are compacted at a lower rate than the InSAR detection threshold accuracy of 1~4 mm/year. Third, the measuring principle of groundwater depl GRACE and the official groundwater budgets is different. GRACE provides water changes consisting of all dynamic changes in aquifers, but it is contaminated by errors because of the inherent resolution of the gravity field. Fourth, the official water budget ignores the dynamic impacts of aquifer depletion, wastewater recha water infiltration. It is worth noting that InSAR-derived displacements contain versible and mainly irreversible. Mass changes in phreatic aquifers perform greate on the compressibility of soil than those of confined aquifers [73]. The groundwa data from unconfined aquifers were utilized in this paper, despite the fact that them contribute to the subsidence. The diversity between simulated GWSA and in situ monitoring wells can be explained from four aspects. At first, GRACE-GWSA does not fully reflect the regional groundwater depletion due to leakage errors. Second, 36% of grids are not calculated in the GWSA trends map at 0.25 • original resolution ((406-258)/406 ≈ 0.36), which means that several grids are compacted at a lower rate than the InSAR detection threshold with an accuracy of 1~4 mm/year. Third, the measuring principle of groundwater depletion by GRACE and the official groundwater budgets is different. GRACE provides water storage changes consisting of all dynamic changes in aquifers, but it is contaminated by leakage errors because of the inherent resolution of the gravity field. Fourth, the official groundwater budget ignores the dynamic impacts of aquifer depletion, wastewater recharge, and water infiltration. It is worth noting that InSAR-derived displacements contain both reversible and mainly irreversible. Mass changes in phreatic aquifers perform greater effects on the compressibility of soil than those of confined aquifers [73]. The groundwater level data from unconfined aquifers were utilized in this paper, despite the fact that both of them contribute to the subsidence.

Conclusions
Both InSAR and GRACE are complementary and sensitivities to aquifer system processes. InSAR technology can be applied to downscale the groundwater storage anomalies retrieved by GRACE in space. The NSFWC scheme is conducive to generating highresolution GWSA from 0.5 • to 0.05 • , indicating the synergy of two remote sensing tools to survey the sustainable development of groundwater. In this paper, the datasets of GWSA (from GRACE and CLSM), land subsidence (from InSAR and GPS), and in situ monitoring wells are jointly analyzed, and we obtained the consistent distribution of groundwater changes in the Beijing Plain. InSAR results manifest that there are currently four cones of depression, and other regions are basically stable during 2017-2020. The driving factors of subsidence are complex and strongly related to the natural consolidation of soil and seasonal changes in precipitation. In general, the long-term component is related to groundwater exploitation, and the seasonal component is affected by climate change, such as precipitation and temperature. On the contrary, it is beneficial to recover the groundwater level via the increased recharge of the SNWD project and the decreased groundwater extraction. At present, it is essential for water authorities to focus on monitoring the regions that are still being overdrawn.
Furthermore, sedimentation grids are filtered from InSAR deformation, and the high correlation between vertical displacements and groundwater level change is verified by their trend components, despite the hydrogeology and geodesy properties. Therefore, InSAR-based deformations are assigned to GRACE-based GWSA as weight coefficients so that the simulated GWSA can offer high-resolution signals for intense groundwater usage. The NSFWC scheme offers a unique perspective on the groundwater issue. The surface displacement detected by satellites is a useful indicator for surveying GWSA. The effectiveness of the simulated method has been proven to enhance the spatial resolution of GWSA estimated by GRACE.
This work is the first to merge the InSAR time-series vertical deformation into GRACEderived GWSA over the Beijing Plain. The simulated high-resolution GWSA data can provide guidance for groundwater management to the sustainable development of water resources. More work toward processing the land subsidence over a longer period of time is needed in the future. Then, spatial downscaling of GRACE combined long-term InSAR data needs to be further explored using machine learning methods (e.g., support vector machine). Additionally, it is interesting to reconstruct the complete GRACE time-series data integrating InSAR data to fill the gap between the GRACE and GRACE-FO missions.