Next Article in Journal
Tracing COVID-19 Trails in Wastewater: A Systematic Review of SARS-CoV-2 Surveillance with Viral Variants
Next Article in Special Issue
The Influence of the South-to-North Water-Diversion Project on Terrestrial Water-Storage Changes in Hebei Province
Previous Article in Journal
Estimation of Lake Storage Based on the Surrounding Topography around the Lake from the SRTM DEM
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
School of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 211100, China
2
China Academy of Aerospace Science and Technology Innovation, Beijing 100076, China
3
Ministry of Ecology and Environment Center for Satellite Application on Ecology and Environment, Beijing 100089, China
4
School of Earth Sciences and Engineering, Hohai University, Nanjing 211100, China
5
School of Geomatics, Liaoning Technical University, Fuxin 123000, China
*
Author to whom correspondence should be addressed.
These authors contributed equally to this paper.
Water 2023, 15(6), 1017; https://doi.org/10.3390/w15061017
Submission received: 2 February 2023 / Revised: 26 February 2023 / Accepted: 6 March 2023 / Published: 7 March 2023
(This article belongs to the Special Issue Application of GRACE Observations in Water Cycle and Climate Change)

Abstract

:
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 reflect 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 coefficient 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.

1. 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 time-consuming and labor-intensive. Thus, drilling observations are insufficient for capturing the spatial details of groundwater variability. A groundwater budget is a common estimation method used by management organizations, but it ignores both aquifer dynamics and local water circulation, which change with GWSA. The shortage of climatic agreements and the difficulty of calculating water fluxes (recharge and discharge) also make it hard to compare different water administrations. Hence, efficient survey tools are necessary to measure the growing pressure on groundwater resources. Recently, remote sensing monitoring technologies such as GRACE and InSAR have shown powerful potential in hydrological applications, which can achieve continuous real-time monitoring.
The Gravity Recovery and Climate Experiment (GRACE) satellite mission was launched in 2002 to detect the mass loss caused by water migration. However, after 15 years (2002–2017) of operation in orbit, GRACE no longer provided available measurement data. Subsequently, the GRACE Follow-On (GRACE-FO) mission was launched by NASA and the German Aerospace Centre (GFZ) to produce more accurate data from June 2018 to now. There is a one-year data gap between GRACE and GRACE-FO, and many previous researchers reconstructed the GRACE time series to fill in the gap [8]. GRACE mascon solutions were released by the Jet Propulsion Laboratory (JPL), the University of Texas at Austin Center for Space Research (CSR), and the NASA Goddard Space Flight Center (GSFC). At present, the GRACE mascon has been widely applied for evaluating hydrogeological and climatological issues. For instance, Alshehri and Mohamed (2023) used GRACE data with the results of the Global Land Data Assimilation System (GLDAS) model to estimate changes in the Wadi As-Sirhan Basin’s groundwater storage and the controlling factors [9]. Mohamed et al. (2023) integrated three GRACE mascon products with GLDAS models to estimate the spatiotemporal changes of TWSA, GWSA, the recharge rates resulting from climate change, and human activities in the Republic of Chad [10].
InSAR technology offers a novel viewpoint to the global groundwater monitoring system, which is increasingly used to monitor ground subsidence and infer hydrogeological properties [11,12]. Previous studies reported that InSAR-based deformation could deduce elastic aquifer storage coefficient, estimate groundwater level (GWL) changes, and calibrate groundwater numerical models [13,14,15]. InSAR detects land deformation because of the soil-layer response to pressure (water head) variations. This method has the advantages of global extent, high precision, and high resolution, but there are certain limitations in regions covered by abundant vegetation. Additionally, the compaction thickness varies in space because of the uneven aquifer response to hydrostatic pressure changes. Even if stress changes seem to impact the structure of compressed layers, some depleted layers are not affected by land displacement greater than the typically detected values of InSAR. Its applications in the evaluation of water resources are constricted to aquifers, which are susceptible to evident motion. Persistent Scatterers Interferometry (PS-InSAR) [16] and Small Baseline Subset (SBAS) [17] algorithms have been widely utilized in urban deformation monitoring [18]. The phase and amplitude of persistent scatterers are relatively stable for a long time, which can effectively reduce the impacts of spatiotemporal de-coherence, atmospheric delay, and noise. The increasing availability of Sentinel-1 images makes it possible to obtain land subsidence data at the basin scale.
GRACE data are sensitive to water resource changes. It offers an effective approach to derive GWSA, supposing that other water components have been estimated and reduced [19,20,21]. However, the footprint of GRACE is approximately 300 km [22], which is applicable to a basin above 100,000 km2. The original resolution of GRACE is inadequate to favor estimation on a regional scale. Many attempts were conducted to enhance GRACE resolution, focusing on new algorithm development [23,24]. The improved GRACE inversion algorithm is significantly better for evaluating groundwater changes by considering the spatial heterogeneity of masses [25,26]. Thus, combining groundwater-related GRACE signals with InSAR measurements can tender high-resolution groundwater depletion data.
The coarse spatial resolution of GRACE has motivated a number of studies focused on downscaling GRACE measurements by both dynamic and statistical methods for use at smaller spatial scales [27,28]. Dynamic downscaling methods downscale GRACE by assimilating 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 InSAR 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.

2. Materials and Methods

2.1. 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.
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.

2.2. Materials

2.2.1. 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.

2.2.2. 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.

2.2.3. 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].

2.3. 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.

2.3.1. 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].
ϕ d e f o = ϕ i n t ϕ r e f + ϕ d e m + ϕ a t m + ϕ r e s , ϕ i n t = 4 π λ Δ R , ϕ r e f = 4 π λ B 0 , ϕ d e m = 4 π λ B 0 h R sin θ ,
where ϕ i n t refers to the interference phase; ϕ r e f stands for the reference ellipsoid phase; ϕ d e m suggests the topographical phase, which refers to the relative ground motion during the two shooting dates. ϕ a t m is the atmospheric delay phase, and ϕ r e s 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 multi-look processing of 10 × 2 was used to improve the signal/noise ratio. Secondly, DEM data were employed to reduce the ϕ d e m 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].
D v x , y = D L O S x , y + sin θ cos φ D n x , y + sin θ sin φ D e x , y δ L O S cos θ , D n x , y = D L O S / cos α , D e x , y = D L O S / sin α
where DLOS is the LOS deformation; Dv is the perpendicular movement; Dn and De indicate the horizontal movements; θ , α , and φ represent radar incidence angle, heading angle, and azimuth angle, respectively; δ L O S is radar noise.

2.3.2. 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 Moisture Storage Anomalies (SMSA) and Surface Water Storage Anomalies (SWSA) from TWSA under the same resolution, which is expressed by the following equation [52,53].
G W S A k = T W S A k S M S A k S W S A k ,
where TWSAk is the terrestrial water storage anomalies of the kth grid from GRACE; SMSAk is the NOAH-derived soil moisture storage anomalies of the kth grid; SWSAk 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].

2.3.3. 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.
First, the deformation estimated by InSAR was re-sampled from vector pixels (50 × 40 m) to raster grids (0.05° × 0.05°) using the kriging spatial interpolation method [56].
Z λ 0 , φ 0 , h 0 = i = 1 n σ i Z λ i , φ i , h i ,
The pending weight coefficient was obtained by the corresponding minimum kriging variance on the premise of unbiased, which was determined by the following equation.
C λ 0 , φ 0 , h 0 , λ j , φ j , h j = i = 1 n σ i C λ i , φ i , h i , λ j , φ j , h j + μ i = 1 n σ i = 1 ,
where σ i stands for the pending weight coefficient; Z λ 0 , φ 0 , h 0 means deformation results based on the kriging interpolation method; Z λ i , φ i , h i is deformation at the coordinate of λ i , φ i , h i , C λ 0 , φ 0 , h 0 , λ j , φ j , h j represents the covariance of Z λ i , φ i , h i and Z λ j , φ j , h j .
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].
D s u b = F i l t e r D v < 0 , N ,
where Dsub is the rate value of decreasing grid; Dv 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]:
y = x x min x max x min ,
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.
W k = 1 D s u b k D v A + D v B + D v C + D v D / n D s u b max D v A + D v B + D v C + D v D / n ,
where Wk is the normalization result of subsidence grids; D v A , D v B , D v C , and D v D are the vertical displacements at four cones of depression (A, B, C, and D) respectively; D s u b k denotes the subsidence rate of the kth grid; D s u b max 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.
G W S A M k = G W S A k × W k + G W S A k ,
where G W S A M k represents the GWSA estimated by the feature combined scheme; GWSAk 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].
G W S A O k = Δ G W L k × S y ,
where G W S A O k suggests the GWSA estimated by groundwater level changes of the kth grid, Δ G W 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 (ek) [57], Pearson correlation coefficient (r) [61], and root mean square error (RMSE) [62]:
e k = G W S A M k G W S A O k ,
r = k = 1 m G W S A M k G W S A M ¯ G W S A O k G W S A O ¯ k = 1 m G W S A M k G W S A M ¯ 2 G W S A O k G W S A O ¯ 2 ,
R M S E = k = 1 m ( G W S A M k G W S A O k ) 2 m ,
where ek denotes the absolute error at the kth grid, G W S A M k indicates the simulation result at the kth grid, G W S A O k means the in situ data at the kth grid, G W S A M ¯ is the average of simulated results, G W S A O ¯ is the average of in situ data. m is the number of settlement grids.

3. Results

3.1. Groundwater Level Changes

Figure 3 and Figure 4 illustrate the spatial distribution of Δ G W L trends and the average sequence of groundwater depth changes in the Beijing Plain, respectively. Figure 3a–c show the Δ G W L trend during 2005–2011, 2012–2014, and 2015–2018, respectively. GWL in 2005–2011 mostly shows decreased trends due to continuous groundwater overexploitation. Then, it becomes gradually stable owing to the increase in rainfall during 2012–2014. Since 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 Δ G W 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.

3.2. InSAR-Derived Land Deformation

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.

3.3. Reliability of InSAR Results

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.

3.4. GWSA Derived by NSFWC Scheme

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.
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.
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.

4. Discussion

4.1. Comparison of GWSA, Land Subsidence, and Precipitation

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.
Figure 12b compares the average subsidence values of InSAR in northern Beijing Plain with GPS at the site of BJSH. The long-term trends and subsidence magnitude of InSAR and GPS are consistent with each other, which suggests that InSAR can capture high-precision deformation. 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.
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.

4.2. 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.

4.3. 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.
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.

5. 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 high-resolution 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 GRACE-derived 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.

Author Contributions

Conceptualization and methodology, Q.W., W.Z. and W.Y.; writing—original draft preparation, Q.W.; writing—review and editing, W.Y., W.Z., G.K., Q.H. and Y.S.; supervision, W.Z., W.Y. and G.K.; funding acquisition, W.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded in part by the National Natural Science Foundation of China, grant number 42274119, in part by the Liaoning Revitalization Talents Program, grant number XLYC2002082, in part by National Key Research and Development Plan Key Special Projects of Science and Technology Military Civil Integration, grant number 2022YFF1400500.

Data Availability Statement

The JPL mascon data are available at https://grace.jpl.nasa.gov/data/get-data/jpl_global_mascons/#opennewwindow (accessed on 17 May 2022). The CSR mascon data are available at http://www.csr.utexas.edu/grace (accessed on 5 May 2022). The NOAH model is available at https://disc.gsfc.nasa.gov/datasets?keywords=GLDAS&page=1 (accessed on 2 May 2022). The Sentinel-1A images are available at https://scihub.copernicus.eu (accessed on 11 March 2021).

Acknowledgments

The editor and reviewers are sincerely acknowledged for their instructive and detailed comments on this manuscript. The authors appreciate all centers for supplying open-source data sets used in the present research. Qingqing Wang, Wei Zheng, and Wenjie Yin contributed equally to this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhang, J.; Sun, Q.; Fourie, A.; Ju, F.; Dong, X. Risk assessment and prevention of surface subsidence in deep multiple coal seam mining under dense above-ground buildings: Case study. Hum. Ecol. Risk Assess. 2019, 25, 1579–1593. [Google Scholar] [CrossRef]
  2. Smith, R.G.; Knight, R.; Chen, J.; Reeves, J.A.; Zebker, H.A.; Farr, T.; Liu, Z. Estimating the permanent loss of groundwater storage in the southern San Joaquin Valley, California. Water Resour. Res. 2017, 53, 2133–2148. [Google Scholar] [CrossRef]
  3. Zhang, Z.; Wang, C.; Wang, M.; Wang, Z.; Zhang, H. Surface deformation monitoring in Zhengzhou city from 2014 to 2016 using time-series InSAR. Remote Sens. 2018, 10, 1731. [Google Scholar] [CrossRef] [Green Version]
  4. Chen, B.; Gong, H.; Lei, K.; Li, J.; Zhou, C.; Gao, M.; Guan, H.; Lv, W. Land subsidence lagging quantification in the main exploration aquifer layers in Beijing plain, China. Int. J. Appl. Earth Obs. 2019, 75, 54–67. [Google Scholar] [CrossRef]
  5. Zhu, Z.; Ren, C.; Zhou, L.; Shi, X.; Li, X.; Zhang, D. Monitoring Tianjin land subsidence by SBAS-InSAR based on Sentinel-1A SAR images. ISPRS J. Photogramm. 2020, 3, 315–331. [Google Scholar] [CrossRef] [Green Version]
  6. Hao, R.; Wang, X.W.; Xu, Y.S. Analysis of the factors resulting in the acceleration of land subsidence in the central area of Dezhou city, China. Environ. Earth Sci. 2022, 81, 400. [Google Scholar] [CrossRef]
  7. Bai, L.; Jiang, L.; Zhao, Y.; Li, Z.; Cao, G.; Zhao, C.; Liu, R.; Wang, H. Quantifying the influence of long-term overexploitation on deep groundwater resources across Cangzhou in the North China Plain using InSAR measurements. J. Hydrol. 2022, 605, 127368. [Google Scholar] [CrossRef]
  8. Liu, B.; Zou, X.; Yi, S.; Sneeuw, N.; Li, J.; Cai, J. Reconstructing GRACE-like time series of high mountain glacier mass anomalies. Remote Sens. Environ. 2022, 280, 113177. [Google Scholar] [CrossRef]
  9. Alshehri, F.; Mohamed, A. Analysis of Groundwater Storage Fluctuations Using GRACE and Remote Sensing Data in Wadi As-Sirhan, Northern Saudi Arabia. Water 2023, 15, 282. [Google Scholar] [CrossRef]
  10. Mohamed, A.; Abdelrady, A.; Alarifi, S.S.; Othman, A. Geophysical and Remote Sensing Assessment of Chad’s Groundwater Resources. Remote Sens. 2023, 15, 560. [Google Scholar] [CrossRef]
  11. Chaussard, E.; Burgmann, R.; Shirzaei, M.; Fielding, E.J.; Baker, B. Predictability of hydraulic head changes and characterization of aquifer system and fault properties from InSAR-derived ground deformation. J. Geophys. Res. Solid Earth 2014, 119, 6572–6590. [Google Scholar] [CrossRef]
  12. Ojha, C.; Shirzaei, M.; Werth, S.; Argus, D.F.; Farr, T.G. Sustained groundwater loss in California’s Central Valley exacerbated by intense drought periods. Water Resour. Res. 2018, 54, 4449–4460. [Google Scholar] [CrossRef]
  13. Castellazzi, P.; Martel, R.; Galloway, D.L. Assessing groundwater depletion and dynamics using GRACE and InSAR: Potential and limitations. Groundwater 2016, 54, 768. [Google Scholar] [CrossRef] [Green Version]
  14. Béjar-Pizarro, M.; Ezquerro, P.; Herrera, G.; Tomás, R.; Guardiola-Albert, C.; Hernández, J.; Merodo, J.; Marchamalo, M.; Martínez, M. Mapping groundwater level and aquifer storage variations from InSAR measurements in the Madrid aquifer, Central Spain. J. Hydrol. 2017, 547, 678–689. [Google Scholar] [CrossRef] [Green Version]
  15. Jiang, L.; Bai, L.; Zhao, Y.; Cao, G.; Wang, H.; Sun, Q. Combining InSAR and hydraulic head measurements to estimate aquifer parameters and storage variations of confined aquifer system in Cangzhou, North China Plain. Water Resour. Res. 2018, 54, 8234–8252. [Google Scholar] [CrossRef]
  16. Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans. Geosci. Remote Sens. 2001, 39, 8–20. [Google Scholar] [CrossRef]
  17. Berardino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375–2383. [Google Scholar] [CrossRef] [Green Version]
  18. Wang, Y.; Huang, Q.; Guo, J.; Jiang, J.; Wang, Q.; Nishyirimbere, N. Displacement monitoring and modeling of Qinghai–Tibet Railway in permafrost area using Sentinel-1A data. Arab. J. Geosci. 2019, 12, 204. [Google Scholar] [CrossRef]
  19. Scanlon, B.R.; Longuevergne, L.; Long, D. Ground referencing GRACE satellite estimates of groundwater storage changes in the California Central Valley, USA. Water Resour. Res. 2012, 48, W04520. [Google Scholar] [CrossRef] [Green Version]
  20. Zhong, Y.; Zhong, M.; Feng, W.; Zhang, Z.; Shen, Y.; Wu, D. Groundwater depletion in the west Liaohe river basin, china and its implications revealed by GRACE and in situ measurements. Remote Sens. 2018, 10, 493. [Google Scholar] [CrossRef] [Green Version]
  21. Yin, W.; Hu, L.; Zheng, W.; Jiao, J.; Han, S.; Zhang, M. Assessing underground water exchange between regions using GRACE data. J. Geophys. Res. Atmos. 2020, 125, e2020JD032570. [Google Scholar] [CrossRef]
  22. Farinotti, D.; Longuevergne, L.; Moholdt, G.; Duethmann, D.; Molg, T.; Bolch, T.; Vorogushyn, S.; Guntner, A. Substantial glacier mass loss in the Tien Shan over the past 50 years. Nat. Geosci. 2015, 8, 716–722. [Google Scholar] [CrossRef]
  23. Yin, W.; Hu, L.; Zhang, M.; Wang, J.; Han, S. Statistical downscaling of GRACE-derived groundwater storage using ET data in the North China Plain. J. Geophys. Res. Atmos. 2018, 123, 5973–5987. [Google Scholar] [CrossRef]
  24. Vishwakarma, B.; Zhang, J.; Sneeuw, N. Downscaling GRACE total water storage change using partial least squares regression. Sci. Data 2021, 8, 95. [Google Scholar] [CrossRef]
  25. Long, D.; Chen, X.; Scanlon, B.R.; Wada, Y.; Hong, Y.; Singh, V.P.; Chen, Y.; Wang, C.; Han, Z.; Yang, W. Have GRACE satellites overestimated groundwater depletion in the Northwest India Aquifer? Sci. Rep. 2016, 6, 24398. [Google Scholar] [CrossRef] [Green Version]
  26. Liu, Z.; Liu, P.; Massoud, E.; Farr, T.; Lundgren, P.; Famiglietti, J. Monitoring Groundwater change in California’s central valley using Sentinel-1 and GRACE observations. Geosciences 2019, 9, 436. [Google Scholar] [CrossRef] [Green Version]
  27. Milewski, A.; Thomas, M.; Seyoum, W.; Rasmussen, T. Spatial downscaling of grace TWSA data to identify spatiotemporal groundwater level trends in the upper Floridian aquifer, Georgia, USA. Remote Sens. 2019, 11, 2756. [Google Scholar] [CrossRef] [Green Version]
  28. Chen, Z.; Zheng, W.; Yin, W.; Li, X.; Zhang, G.; Zhang, J. Improving the Spatial Resolution of GRACE-Derived Terrestrial Water Storage Changes in Small Areas Using the Machine Learning Spatial Downscaling Method. Remote Sens. 2021, 13, 4760. [Google Scholar] [CrossRef]
  29. Zhong, D.; Wang, S.; Li, J. A self-calibration variance-component model for spatial downscaling of grace observations using land surface model outputs. Water Resour. Res. 2020, 11, 376–382. [Google Scholar] [CrossRef]
  30. Arshad, A.; Zhang, W.; Zhang, Z.; Wang, S.; Shalamzari, M.J. Reconstructing high-resolution gridded precipitation data using an improved downscaling approach over the high altitude mountain regions of Upper Indus Basin (UIB). Sci. Total Environ. 2021, 784, 147140. [Google Scholar] [CrossRef]
  31. Arshad, A.; Ali, M.; Maryam, S.; Ahmad, B. Combining downscaled-GRACE data with SWAT to improve the estimation of groundwater storage and depletion variations in the Irrigated Indus Basin (IIB). Sci. Total Environ. 2022, 838, 156044. [Google Scholar] [CrossRef]
  32. Ali, S.; Liu, D.; Fu, Q.; Cheema, M.; Pal, S.; Arshad, A.; Pham, Q.; Zhang, L. Constructing high-resolution groundwater drought at spatio-temporal scale using GRACE satellite data based on machine learning in the Indus Basin. J. Hydrol. 2022, 612, 128295. [Google Scholar] [CrossRef]
  33. Ali, S.; Khorrami, B.; Jehanzaib, M.; Tariq, A.; Ajmal, M.; Arshad, A.; Shafeeque, M.; Dilawar, A.; Basit, I.; Zhang, L.; et al. Spatial Downscaling of GRACE Data Based on XGBoost Model for Improved Understanding of Hydrological Droughts in the Indus Basin Irrigation System (IBIS). Remote Sens. 2023, 15, 873. [Google Scholar] [CrossRef]
  34. Noor, R.; Arshad, A.; Shafeeque, M.; Liu, J.; Baig, A.; Ali, S.; Maqsood, A.; Pham, Q.B.; Dilawar, A.; Khan, S.N.; et al. Combining APHRODITE Rain Gauges-Based Precipitation with Downscaled-TRMM Data to Translate High-Resolution Precipitation Estimates in the Indus Basin. Remote Sens. 2023, 15, 318. [Google Scholar] [CrossRef]
  35. Pranjal, P.; Kadiyan, N.; Chatterjee, R.S. Interpreting land subsidence impacts due to groundwater depletion using remote sensing-based GRACE gravity anomaly and D-InSAR technique: A study on north-western parts of India. Environ. Earth Sci. 2021, 80, 596. [Google Scholar] [CrossRef]
  36. Castellazzi, P.; Longuevergne, L.; Martel, M.; Rivera, A.; Brouard, C.; Chaussard, E. Quantitative mapping of groundwater depletion at the water management scale using a combined GRACE/InSAR approach. Remote Sens. Environ. 2018, 205, 408–418. [Google Scholar] [CrossRef]
  37. Shang, Q.; Liu, X.; Deng, X.; Zhang, B. Downscaling of GRACE datasets based on relevance vector machine using InSAR time series to generate maps of groundwater storage changes at local scale. J. Appl. Rem. Sens. 2019, 13, 048503. [Google Scholar] [CrossRef]
  38. Massoud, E.C.; Liu, Z.; Shaban, A.; Hage, M.E. Groundwater depletion signals in the Beqaa Plain, Lebanon: Evidence from GRACE and Sentinel-1 data. Remote Sens. 2021, 13, 915. [Google Scholar] [CrossRef]
  39. Castellazzi, P.; Garfias, J.; Martel, R.; Brouard, C.; Rivera, A. InSAR to support sustainable urbanization over compacting aquifers: The case of Toluca Valley, Mexico. Int. J. Appl. Earth Obs. Geoinf. 2017, 63, 33–44. [Google Scholar] [CrossRef]
  40. He, Q. Temporal InSAR Monitoring of Surface Deformation and Its Relationship with Groundwater Level in Beijing. Master’s Thesis, China University of Geosciences, Beijing, China, 2020. (In Chinese). [Google Scholar]
  41. Liu, Y.; Ye, C.; Jia, S. Division of water-bearing rock group and compressible layer in subsidence area of Beijing Plain. Urban Geol. 2007, 4, 10–15. (In Chinese) [Google Scholar]
  42. Torres, R.; Snoeij, P.; Geudtner, D.; Bibby, D.; Davidson, M.; Attema, E.; Potin, P.; Rommen, B.; Floury, N.; Brown, M.; et al. GMES Sentinel-1 mission. Remote Sens. Environ. 2012, 120, 9–24. [Google Scholar] [CrossRef]
  43. Wiese, D.N.; Landerer, F.W.; Watkins, M.M. Quantifying and reducing leakage errors in the JPL RL05M GRACE mascon solution. Water Resour. Res. 2016, 52, 7490–7502. [Google Scholar] [CrossRef]
  44. Save, H.; Bettadpur, S.; Tapley, B.D. High-resolution CSR GRACE RL05 MASCONs. J. Geophys. Res. Solid. Earth. 2016, 121, 7547–7569. [Google Scholar] [CrossRef]
  45. Shen, Y.; Zheng, W.; Yin, W.; Xu, A.; Zhu, H.; Yang, S.; Su, K. Inverted Algorithm of Terrestrial Water-Storage Anomalies Based on Machine Learning Combined with Load Model and Its Application in Southwest China. Remote Sens. 2021, 13, 3358. [Google Scholar] [CrossRef]
  46. Crosetto, M.; Monserrat, O.; Cuevas-Gonz´alez, M.; Devanth´ery, N.; Crippa, B. Persistent scatterer interferometry: A review. ISPRS J. Photogramm. 2016, 115, 78–89. [Google Scholar] [CrossRef] [Green Version]
  47. Huang, Q.H.; Monserrat, O.; Crosetto, M.; Crippa, B.; Wang, Y.A.; Jiang, J.F.; Ding, Y. Displacement monitoring and health evaluation of two bridges using Sentinel-1 SAR images. Remote Sens. 2018, 10, 1714. [Google Scholar] [CrossRef] [Green Version]
  48. He, X.; He, M. InSAR Earth Observation Data Processing Method and Comprehensive Survey; Science Press: Beijing, China, 2012; pp. 54–56. [Google Scholar]
  49. Wang, Q.Q.; Huang, Q.H.; He, N.; He, B.; Wang, Z.C.; Wang, Y.A. Displacement monitoring of Upper Atbara Dam based on time series InSAR. Surv. Rev. 2019, 52, 485–496. [Google Scholar] [CrossRef]
  50. Han, P. Status quo of groundwater development and utilization in Haihe River basin and its management. Haihe Water Resour. 2015, 39, 1–5. (In Chinese) [Google Scholar]
  51. Ren, Y.; Pan, Y.; Gong, H. Haihe Basin groundwater reserves space trend analysis. J. Cap. Norm. Univ. 2014, 35, 89–98. (In Chinese) [Google Scholar]
  52. Zhong, Y.; Feng, W.; Humphrey, V.; Zhong, M. Human-induced and climate-driven contributions to water storage variations in the haihe river basin, china. Remote Sens. 2019, 11, 3050. [Google Scholar] [CrossRef] [Green Version]
  53. Mohamad, N.; Ahmad, A.; Din, A. Monitoring groundwater depletion due to drought using satellite gravimetry: A review. IOP Conf. Ser. Earth Environ. Sci. 2020, 540, 012054. [Google Scholar] [CrossRef]
  54. Yang, X.; Tian, S.; Feng, W.; Ran, J.; You, W.; Jiang, Z.; Gong, X. Spatiotemporal evaluation of water storage trends from hydrological models over Australia using GRACE Mascon solutions. Remote Sens. 2020, 12, 3578. [Google Scholar] [CrossRef]
  55. Vautard, R.; Ghil, M. Singular spectrum analysis in nonlinear dynamics with applications paleoclimate time series. Physica D 1989, 35, 399–424. [Google Scholar] [CrossRef]
  56. Li, Y.; Zuo, X.; Xiong, P.; You, H.; Zhang, H.; Yang, F.; Zhao, Y.; Yang, Y.; Liu, Y. Deformation monitoring and analysis of Kunyang phosphate mine fusion with InSAR and GPS measurements. Adv. Space Res. 2022, 69, 2637–2658. [Google Scholar] [CrossRef]
  57. Li, Q. Numerical Analysis; Tsinghua University Press: Beijing, China, 2008; pp. 34–35. (In Chinese) [Google Scholar]
  58. Ioffe, S.; Szegedy, C. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Proceedings of the 32nd International Conference on Machine Learning, Lille, France, 6–11 July 2015; Volume 37, pp. 448–456. [Google Scholar]
  59. Zheng, L.Q.; Pan, Y.; Gong, H.L. Comparing groundwater storage changes in two main grain producing areas in china: Implications for sustainable agricultural water resources management. Remote Sens. 2020, 12, 21–43. [Google Scholar] [CrossRef]
  60. Mohamed, A.; Gonçalvès, J. Hydro-geophysical monitoring of the North Western Sahara Aquifer System’s groundwater resources using gravity data. J. Afr. Earth Sci. 2021, 178, 104188. [Google Scholar] [CrossRef]
  61. Rodgers, J.L.; Nicewander, W.A. Thirteen ways to look at the correlation coefficient. Am. Stat. 1988, 42, 59–66. [Google Scholar] [CrossRef]
  62. Chai, T.; Draxler, R.R. Root mean square error (RMSE) or mean absolute error (MAE)?—Arguments against avoiding RMSE in the literature. Geosci. Model. Dev. 2014, 7, 1247–1250. [Google Scholar] [CrossRef] [Green Version]
  63. Li, D.; Hu, L.; Zhang, M.; Liu, X. Evaluation of groundwater withdrawals in Tongzhou district, Beijing, China. Hum. Ecol. Risk Assess. 2019, 27, 1441–1457. [Google Scholar] [CrossRef]
  64. Dong, J.; Lai, S.; Wang, N.; Wang, Y.; Zhang, L.; Liao, M. Multi-scale deformation monitoring with Sentinel-1 InSAR analyses along the Middle Route of the South-North Water Diversion Project in China. Int. J. Appl. Earth Obs. 2021, 100, 102324. [Google Scholar] [CrossRef]
  65. Ge, D.Q. Study on Key Techniques of InSAR Monitoring of Regional Land Subsidence. Ph.D. Thesis, China University of Geosciences, Beijing, China, 2013. (In Chinese). [Google Scholar]
  66. Zhang, T.; He, Q.; Jing, C.; Liu, L. Dynamic monitoring of groundwater in Beijing Plain area based on InSAR. Sci. Technol. Eng. 2019, 19, 16–22. (In Chinese) [Google Scholar]
  67. Jiang, Y.; Tian, F.; Luo, Y. Relationship between stratified land subsidence and groundwater level variation in typical areas of Beijing. South--North Water Transf. Water Sci. Technol. 2015, 13, 95–99. (In Chinese) [Google Scholar]
  68. Hoffmann, J.; Galloway, D.L.; Zebker, H.A. Inverse modeling of interbred storage parameters using land subsidence observations, Antelope Valley, California. Water Resour. Res. 2003, 39, 1031. [Google Scholar] [CrossRef]
  69. Hu, X.; Lu, Z.; Wang, T. Characterization of hydrogeological properties in Salt Lake Valley, Utah using InSAR. J. Geophys. Res. 2018, 123, 1257–1271. [Google Scholar] [CrossRef] [Green Version]
  70. Galloway, D.L.; Hoffmann, J. The application of satellite differential SAR interferometry-derived ground displacements in hydrogeology. Hydrgeol. J. 2007, 15, 133–154. [Google Scholar] [CrossRef] [Green Version]
  71. Hanssen, R.F. Radar Interferometry: Data Interpretation and Error Analysis; Springer Science and Business Media: New York, NY, USA, 2001; pp. 120–121. [Google Scholar]
  72. Ghil, M.; Allen, M.; Dettinger, M. Advanced spectral methods for climatic time series. Rev. Geophys. 2002, 40, 1003. [Google Scholar] [CrossRef] [Green Version]
  73. Pool, D.R. The utility of gravity and water-level monitoring at alluvial aquifer wells in southern Arizona. Geophysics 2008, 73, 49–59. [Google Scholar] [CrossRef]
Figure 1. The geographical overview of the study area.
Figure 1. The geographical overview of the study area.
Water 15 01017 g001
Figure 2. The flowchart of data processing. (a) 0.25° CSR-GWSA; (b) 0.05° CSR-GWSA; (c) 0.5° CSR-GWSA; (d) 0.05° simulated CSR-GWSA.
Figure 2. The flowchart of data processing. (a) 0.25° CSR-GWSA; (b) 0.05° CSR-GWSA; (c) 0.5° CSR-GWSA; (d) 0.05° simulated CSR-GWSA.
Water 15 01017 g002
Figure 3. The spatial distribution of △GWL trends during (a) 2005–2011, (b) 2012–2014, and (c) 2015–2018, respectively.
Figure 3. The spatial distribution of △GWL trends during (a) 2005–2011, (b) 2012–2014, and (c) 2015–2018, respectively.
Water 15 01017 g003
Figure 4. Average △GWL in the Beijing Plain from 2005 to 2018 and GWL trends in different periods.
Figure 4. Average △GWL in the Beijing Plain from 2005 to 2018 and GWL trends in different periods.
Water 15 01017 g004
Figure 5. The GWL changes of ten downtrend wells during 2015–2018.
Figure 5. The GWL changes of ten downtrend wells during 2015–2018.
Water 15 01017 g005
Figure 6. Deformation retrieved by InSAR in the Beijing Plain from May 2017 to July 2020. (a) Vector image with persistent scatterers. The cones of depression are located at A, B, C, and D, respectively; (b) Raster image at the resolution of 0.05°.
Figure 6. Deformation retrieved by InSAR in the Beijing Plain from May 2017 to July 2020. (a) Vector image with persistent scatterers. The cones of depression are located at A, B, C, and D, respectively; (b) Raster image at the resolution of 0.05°.
Water 15 01017 g006
Figure 7. Subsidence sequences at four funnels from May 2017 to July 2020. (a) Cumulative subsidence; (b) Average subsidence rate.
Figure 7. Subsidence sequences at four funnels from May 2017 to July 2020. (a) Cumulative subsidence; (b) Average subsidence rate.
Water 15 01017 g007
Figure 8. Sedimentation grids extracted from InSAR deformation.
Figure 8. Sedimentation grids extracted from InSAR deformation.
Water 15 01017 g008
Figure 9. Comparison among InSAR, GRACE, CLSM, and △GWL sequence during 2017–2020. (ad) The original sequence comparison at the in situ points of P17, P66, P68, and P86, respectively; (eh) The correlations of trend components between subsidence and groundwater changes at the four in situ points.
Figure 9. Comparison among InSAR, GRACE, CLSM, and △GWL sequence during 2017–2020. (ad) The original sequence comparison at the in situ points of P17, P66, P68, and P86, respectively; (eh) The correlations of trend components between subsidence and groundwater changes at the four in situ points.
Water 15 01017 g009
Figure 10. Annual trends of variables retrieved with GRACE and GLDAS during 2017–2020. (a) 0.5° JPL-TWSA; (b) 0.25° CSR-TWSA; (c) 0.25° NOAH-SMSA; (d) 0.25° NOAH-SWSA; (eh) the interpolated results at 0.05° corresponding to (ad), respectively. The grey dots represent areas with non-significant trends (p ≥ 0.05).
Figure 10. Annual trends of variables retrieved with GRACE and GLDAS during 2017–2020. (a) 0.5° JPL-TWSA; (b) 0.25° CSR-TWSA; (c) 0.25° NOAH-SMSA; (d) 0.25° NOAH-SWSA; (eh) the interpolated results at 0.05° corresponding to (ad), respectively. The grey dots represent areas with non-significant trends (p ≥ 0.05).
Water 15 01017 g010
Figure 11. The simulated GWSA based on the NSFWC scheme. (a) 0.5° JPL-GWSA; (b) 0.5° CSR-GWSA; (c) 0.5° CLSM-GWSA; (d) 0.05° simulated JPL-GWSA; (e) 0.05° simulated CSR-GWSA; (f) 0.05° simulated CLSM-GWSA. The grey dots represent areas with non-significant trends (p ≥ 0.05).
Figure 11. The simulated GWSA based on the NSFWC scheme. (a) 0.5° JPL-GWSA; (b) 0.5° CSR-GWSA; (c) 0.5° CLSM-GWSA; (d) 0.05° simulated JPL-GWSA; (e) 0.05° simulated CSR-GWSA; (f) 0.05° simulated CLSM-GWSA. The grey dots represent areas with non-significant trends (p ≥ 0.05).
Water 15 01017 g011
Figure 12. The relationship among GWSA, subsidence, and precipitation. (a) Comparison of average values in the Beijing Plain; (b) Comparison of average values in northern Beijing Plain; (c) Comparison of seasonal components.
Figure 12. The relationship among GWSA, subsidence, and precipitation. (a) Comparison of average values in the Beijing Plain; (b) Comparison of average values in northern Beijing Plain; (c) Comparison of seasonal components.
Water 15 01017 g012
Figure 13. Comparison of subsidence sequence between individual grids nearby site BJSH and average values in the northern Beijing Plain.
Figure 13. Comparison of subsidence sequence between individual grids nearby site BJSH and average values in the northern Beijing Plain.
Water 15 01017 g013
Figure 14. The topography of four depression funnels in Google Earth. (a), (b), (c) and (d) shows the topography of A, B, C, and D, respectively. The location of A, B, C, and D can be found in Figure 6.
Figure 14. The topography of four depression funnels in Google Earth. (a), (b), (c) and (d) shows the topography of A, B, C, and D, respectively. The location of A, B, C, and D can be found in Figure 6.
Water 15 01017 g014
Figure 15. The errors (mm/year) of GWSA trends from the residual between simulated results and △GWL. The region outside the Beijing Plain is set with a transparency of 30% in the color map.
Figure 15. The errors (mm/year) of GWSA trends from the residual between simulated results and △GWL. The region outside the Beijing Plain is set with a transparency of 30% in the color map.
Water 15 01017 g015
Table 1. Parameters of Sentinel-1A sensor and datasets used in this study.
Table 1. Parameters of Sentinel-1A sensor and datasets used in this study.
Sentinel-1AParameterDatasetParameter
Launch time2014.4.3Sensing period20 May 2017–15 July 2020
BandC-band (~5.6 cm)Image number39
Coverage width250 kmOrbit number142
Spatial resolution5 m × 20 mMulti-look10 × 2
Track height700 kmOrderAscending
Revisit time12 daySensor modeIW
Sensor modeSM, IW, EW, WVPolarisationVV
PolarisationVV, HH,
VV+VH, HH+HV
Time baseline24/36 days
Space baseline<180 m
Table 2. Correlation coefficients between subsidence and △GWL at the four in situ points (P17, P66, P68, and P86).
Table 2. Correlation coefficients between subsidence and △GWL at the four in situ points (P17, P66, P68, and P86).
IDP17P66P68P86
r1 (InSAR/△GWL)0.970.980.940.97
r2 (InSAR/CLSM)0.950.970.940.86
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Wang, Q.; Zheng, W.; Yin, W.; Kang, G.; Huang, Q.; Shen, Y. Improving the Resolution of GRACE/InSAR Groundwater Storage Estimations Using a New Subsidence Feature Weighted Combination Scheme. Water 2023, 15, 1017. https://doi.org/10.3390/w15061017

AMA Style

Wang Q, Zheng W, Yin W, Kang G, Huang Q, Shen Y. Improving the Resolution of GRACE/InSAR Groundwater Storage Estimations Using a New Subsidence Feature Weighted Combination Scheme. Water. 2023; 15(6):1017. https://doi.org/10.3390/w15061017

Chicago/Turabian Style

Wang, Qingqing, Wei Zheng, Wenjie Yin, Guohua Kang, Qihuan Huang, and Yifan Shen. 2023. "Improving the Resolution of GRACE/InSAR Groundwater Storage Estimations Using a New Subsidence Feature Weighted Combination Scheme" Water 15, no. 6: 1017. https://doi.org/10.3390/w15061017

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop