Downscaling Simulation of Groundwater Storage in the Beijing, Tianjin, and Hebei Regions of China Based on GRACE Data

Abstract

Gravity Recovery and Climate Experiment (GRACE)-derived groundwater storage anomalies (GWSA) have been used to highlight groundwater depletion in regional aquifer systems worldwide. However, the use of GRACE products in smaller areas is limited owing to the coarse spatial resolution of the data product. This study utilized a dynamic downscaling method to improve the GWSA resolution from 1° to 0.05° by constructing a groundwater storage numerical model in the Beijing, Tianjin, and Hebei regions of China (BTH). The results indicate that: (1) the GRACE-derived and calculated GWSA had a good match with an average root mean squared error (RMSE) of 2.61 cm equivalent water height (EWH) and an average Nash–Sutcliffe efficiency coefficient (NSE) of 0.84 for the calibration period. (2) The hydraulic gradient coefficient and specific yield had the highest sensitivity, and transmissivity had the lowest sensitivity; however, different forcing data had no obvious influence on the GWSA. (3) The downscaled results not only exhibited time series variations that were consistent with those of the GRACE-derived solutions but also revealed a finer spatial heterogeneity of the GWSA along with increasing correlation coefficients between the GRACE-derived GWSA and the in situ measurements of groundwater levels by 0.06 and reducing the RMSE by 8.85%. (4) The downscaled results reflected the spatiotemporal change characteristics of groundwater storage in different hydrogeological units and administrative regions well. This study demonstrates the potential applications of the proposed downscaling method for both regional and local groundwater resource management.

1. Introduction

Groundwater resources are critical for global food and water security [1,2]. In recent decades, excessive groundwater depletion has become a global problem, affecting economic and social development [3,4]. Due to considerable groundwater depletion [5,6], groundwater storage (GWS) changes under the influence of climate change and human activity have garnered increasing attention. Utilizing field and remote sensing observations for the quantitative estimation of the spatiotemporal variations in groundwater storage is critical for the sustainable management of water resources [7].

The Gravity Recovery and Climate Experiment (GRACE) mission and Its successor, GRACE Follow-On (GRACE-FO) have provided a new opportunity to monitor large-scale GWS changes with unprecedented accuracy by vertically disaggregating terrestrial water storage (TWS) change signals since 2002 [8,9,10]. GRACE gravity satellites effectively compensate for the shortcomings of ground-based measurements, including limited monitoring data coverage, uneven spatial distribution, difficulty in data acquisition, and point-specific estimates [11,12]. GRACE data have been used worldwide to estimate groundwater storage change, such as in Northwest India [13], the U.S. High Plains and California’s Central Valley [14,15], the Colorado River Basin [16], Australia [17,18], the Middle East [19], and North China [11,20,21,22,23,24]. In addition, some studies have proved that the GRACE data and hydrological model are mutually independent but complementary tools, that is, the GRACE-derived data can not only be used to evaluate the groundwater storage change at regional scale but can also be used to calibrate the regional groundwater model and reduce the uncertainty of the model [25]. GRACE data [26] are a significant guide for groundwater modeling, and the combined method composed of GRACE observation data, physical modeling and in situ monitoring data can more effectively manage groundwater. Most researchers applied the hydrological model or groundwater flow model [27,28,29] to verify the reliability of GRACE-derived groundwater storage changes. In order to reduce the uncertainty of regional groundwater storage changes estimation and improve the model reliability, GRACE data are also used for parameter calibration and data assimilation [30]. In recent years, with the development of remote sensing technologies, the integration of GRACE data and other multisource datasets has become a trend in water resource analysis [31,32,33,34]. Although remarkable achievements have been made for large areas, the application of GRACE observations is extremely limited in local areas because of the coarse spatial resolution [35]. Therefore, it is necessary to downscale the GRACE data to a higher resolution that could be practical for refined groundwater resource management at local scales.

Statistical and dynamical downscaling are the two most commonly used methods in downscaling applications [36]. Statistical approaches, such as multivariate regression, artificial neural networks, and machine learning, rely on various statistical techniques to approximate the stochastic between large-scale data and data from small-scale observations, which are comparatively cheap and computationally efficient [37]. Ning et al. [38] constructed an empirical regression model based on the water balance equation to improve the spatial resolution of the TWS. Yin et al. [39] used evapotranspiration data with a relatively high spatial resolution to downscale GRACE-derived groundwater storage anomalies (GWSA). Sun [40] downscaled GRACE-derived TWS by developing an artificial neural network model in the United States. Although statistical downscaling has performed well in improving GRACE-derived data resolution, some disadvantages still need to be considered: the rationale behind statistical downscaling models is not yet well understood; these models do not fully account for some fundamental hydrological processes (e.g., antecedent water storage) and the principle of water balance [41]; and they are not applicable to areas where the correlation between large-scale fields and local elements is not obvious [39].

In contrast to statistical downscaling, dynamic downscaling is used to integrate target datasets into numerical models to produce responses based on physically consistent processes using data assimilation methods. Based on specific physical relationships, dynamic downscaling uses global models to drive regional models, which can be applied at smaller scales to produce finer-resolution information without being limited by observations or the degree of correlation between the large-scale and local elements [36]. Simultaneously, the relatively complete physical process and parameterization scheme of the regional model also guarantees the accuracy of the simulation results. For example, Eicker et al. [42] assimilated GRACE-derived TWS into the WaterGAP Global Hydrology Model using a Kalman filter method in the Mississippi River Basin with reasonable results. Zhong et al. [43] downscaled the GRACE-derived TWS using an iterative adjustment method based on the self-calibration variance-component model and improved the TWS uncertainties through a combination of GRACE observations and the TWS simulated using the Land Surface Model. In addition, many studies [44,45,46] have demonstrated that dynamic downscaling methods can downscale the resolution of GRACE products to a model scale. These methods were based mostly on hydrological or land surface models, which regarded the groundwater in a grid as a linear reservoir and thus did not consider the influences of the aquifer parameters.

The North China Plain (NCP) is one of the regions where the relationship between water resources, food, energy, and ecology is the most significant and the joint risk is the most prominent [21]. Many scholars have studied the spatiotemporal changes of GWS in the NCP using such techniques as groundwater flow modelling [47], data analysis [11,20,21,22], combined hydrological models [24,48], and evaluations using monitoring wells [49,50]. Although these studies have explained the spatiotemporal changes of GWS over the region, a detailed evaluation using the GRACE and ground-observation data from municipalities such as Cangzhou, Xingtai, Handan, and Beijing is still not sufficient for the refined and cross-scale management of groundwater resources. Therefore, downscaling research must be conducted to provide high-resolution water storage estimates in the subregions. To the best of our knowledge, there has been a limited focus on groundwater flow models that apply freely available GRACE-derived GWSA data to conduct dynamic downscaling research.

The objective of this study is to construct a groundwater storage model with the freely available GRACE data to downscale the GWSA from 1° to 0.05°. The only difference between the proposed groundwater storage model and the traditional groundwater flow model is that the groundwater storage changes, instead of the hydraulic head, are used as the key variables in the equation, providing a clear physical meaning. In the next section, the study area information and data sources are prepared. In Section 3, the model downscaling method as well as model evaluation methods are described. In Section 4, the model calibration and validation, presentation, and verification of downscaled results are discussed. In Section 5, the downscaling results are used to evaluate the spatiotemporal changes of GWS and compared to the groundwater level (GWL) changes in different subregions; in addition, model parameters and uncertainty are also discussed. The results can be used as guidelines for enhanced groundwater resource management.

2. Materials

2.1. Study Area

The regions of concern in this study are the Beijing, Tianjin, and Hebei (BTH) Regions in China, which includes most of the NCP (Figure 1). The BTH region is the political, cultural, and economic center, extending from 36°03′N to 42°40′N and 113°27′E to 119°50′E, including the capital city Beijing, Tianjin municipality, and 11 city-level administrative entities in Hebei Province, with an area of approximately 218,000 km² [51]. The BTH region is within a littoral and semi-arid climatic zone; nearly 70% of the rainfall is concentrated from June to September. The average annual precipitation of the entire area is approximately 400–600 mm but less than 400 mm in the northwest plateau and hilly areas. The mean annual evapotranspiration of the region ranges from 1100 to 2000 mm.

2.2. Geology and Hydrogeology

The BTH contains a variety of geographical elements, from the Yanshan-Taihang Mountains in the northwest to the southeast, and gradually transitions into the NCP and the coastal plain, with basins and hills interspersed in the middle, exhibiting a terrain decrease from the southeast to the northwest (Figure 1). The plateau is mainly distributed in the northwest region. The BTH Plain, which is backed by mountains and faces the sea, is part of the NCP. According to the relative positions and causes, from west to east, the region can be divided into three principal zones: the piedmont plains (PP), central plains, and coastal plains. The piedmont and central plains have slow water circulation. Some studies have merged the central and coastal plains into the East Central Plain (ECP) [22,50]. The mountainous areas in the western and northwestern areas of the study area are identified as WM in this study.

Quaternary loose sediments with a thickness ranging from 150 to 600 m, are widely distributed throughout the plain. The aquifer structure transits from a single aquifer in the upper part of piedmont in the western fan to a multi-layer aquifer in the east, and the lithology transits from gravel and pebble to sand, silt, and clay. The PP is the primary infiltration and replenishment area for the regional groundwater which has strong permeability. The central plain has a poor ability to receive infiltration from atmospheric precipitation due to contain a thick continuous impermeable layer. The coastal plain has poor aquifer permeability, water storage, and water supply capabilities. According to the burial and hydraulic properties, the groundwater can be divided into four aquifers from top to bottom. The first aquifer is unconfined groundwater, which is mainly represented in the PP with a thickness of approximately 15 to 55 m; the second to fourth layers are confined aquifers, which predominantly comprise the ECP [22]. Groundwater recharge mainly originates from precipitation infiltration, lateral flow in the mountains, and return flow from irrigation, whereas groundwater discharge primarily occurs through groundwater exploitation and phreatic evaporation [52].

The discrepancy between water supply and demand in the BTH region is prominent, and the per capita volume of water resources is less than one ninth of the national average [51]. The extensive shortage of water resources has led to a serious problem of groundwater overexploitation, with the agricultural sector as the main driver in the BTH region. Overexploitation of groundwater is concentrated in the BTH plains, which has resulted in serious ecological deterioration. The South-to-North Water Diversion project was initiated to alleviate this water stress and has been supplying water to the BTH since 2014 [53]. The Action Program for Comprehensive Treatment of Groundwater Overexploitation in the NCP was issued by the China State Council in 2019 and has been adopted in the northern part of the NCP for the comprehensive treatment of groundwater overexploitation.

2.3. Data Sources

We have summarized the data used in this study (

Table 1. Summary of data used in the study.

Type	Product	Spatial (Temporal) Resolution	Period	Source
TWS	GRACE (CSR-RL06M)	0.5° × 0.5° (monthly)	2003–2020	https://grace.jpl.nasa.gov/data/get-data/ (accessed on 11 March 2022).
SMS, SWE	GLDAS	1° × 1° (monthly)	2003–2020	https://disc.gsfc.nasa.gov/datasets (accessed on 11 March 2022).
Precipitation	TRMM 3B43	0.25° × 0.25° (monthly)	2003–2019	https://gpm.nasa.gov/data/sources (accessed on 21 March 2022).
	ERA5	0.25° × 0.25° (monthly)	2003–2020	https://cds.climate.copernicus.eu/cdsapp#!/home (accessed on 23 March 2022).
	CMA	0.5° × 0.5° (monthly)	2003–2020	http://data.cma.cn/ (accessed on 23 March 2022).
	PENG	0.01° × 0.01° (monthly)	2003–2020	http://poles.tpdc.ac.cn/zh-hans/data/faae7605-a0f2-4d18-b28f-5cee413766a2/ (accessed on 2 April 2022).
AET	GLEAM	0.25° × 0.25° (monthly)	2003–2019	https://www.gleam.eu/ (accessed on 16 March 2022).
	ERA5	0.25° × 0.25° (monthly)	2003–2020	https://cds.climate.copernicus.eu/cdsapp#!/home (accessed on 23 March 2022).
	MOD16	0.01° × 0.01° (monthly)	2003–2020	https://search.earthdata.nasa.gov/search?q=MOD16A2+V006 (accessed on 16 March 2022).
Groundwater level	Ground observations	Stations (monthly)	2005–2014	China Geological Environment Monitoring Institute
			2014, 2016, 2018–2020	Ministry of Water Resources

), and detailed descriptions are provided below.

(1)

In situ GWL data

Monthly in situ GWL measurements from 986 monitoring wells were collected from the China Geological Environment Monitoring Institute and Ministry of Water Resources, including shallow (~10 to 55 m below the surface) and deep aquifers (~60 m to several hundred meters below the surface). Different observations have different starting times, and the periods of available data varied from station to station. All the raw well observations were tested for quality during the preprocessing stage. Data with consecutive gaps of more than 2 years were not used. Following this strict quality control, 143 monitoring well stations were selected from the China Geological Environment Monitoring Institute between 2005 to 2014, and 843 monitoring well stations were selected from the Ministry of Water Resources in 2014, 2016, and 2018 to 2020. Note that the GWL measurements were without conversion to groundwater storage due to the uncertainty of the specific yield, which was mainly used for the qualitative validation.

(2)

GRACE-derived GWSA data

In this study, version 2 of the GRACE (January 2003–June 2017) and GRACE-FO (June 2018–December 2020) RL06 mascon solution, processed by the Jet Propulsion Laboratory (JPL), was used to estimate the monthly TWS changes. The JPL mascon solution uses an area concentration function to constrain the mass variation, which is expressed as a 3° × 3° spherical caps of equal global area. Although this solution provides TWS data with a 0.5° × 0.5° spatial sampling resolution, the values of each grid cell in the same cap were identical and were further resampled to 1° × 1° [50]. More information on the data processing of mason solutions and the leakage error reduction of the JPL mascon solution can be found in Watkins et al. [54], Wiese et al. [55], and Save et al. [56]. The GRACE data short-term gaps were filled using the linear difference method during 2003–2020, except for the 11-month gap (July 2017–May 2018) between GRACE and GRACE-FO.

The monthly data of soil moisture storage (SMS) and snow water equivalent (SWE) were obtained from the Global Land Data Assimilation System (GLDAS) Noah V2.1 model with a spatial resolution of 1° × 1° because no obvious differences in SMS and SWE were found in four land surface models (VIC, CLM, CLSM, and Noah) [39]. The SMS and SWE were also estimated by removing the average values between 2004 and 2009 to keep consistency with the GRACE-derived TWS anomalies. Considering the scarcity of surface water in the arid area of the NCP, the surface water storage capacity of rivers, lakes, and reservoirs was neglected [57]. Therefore, we could obtain GWSA by deducting the amount of soil water anomalies and snow anomalies in the area from the TWS anomalies.

(3)

Precipitation and actual evapotranspiration data

The precipitation data were obtained from the Tropical Rainfall Measuring Mission (TRMM) 3B43, China Meteorological Administration (CMA), European Centre for Medium-Range Weather Forecasts ERA5 reanalysis, and PENG datasets [58]. Actual evapotranspiration (AET) datasets could be obtained from the Global Land-surface Evaporation: the Amsterdam Methodology (GLEAM) V3.5 dataset, ERA5 reanalysis, and Moderate-Resolution Imaging Spectroradiometer (MODIS) ET algorithm (MOD16) AET datasets. The different datasets are compared in Figure 2. In the precipitation products, the TRMM and ERA5 data were similar, while the PENG datasets were relatively small. The AET products varied greatly in summer, and MOD16 overestimated the values from May to August but underestimated the values in other months. The correlations between the different rainfall and evapotranspiration datasets were all close to 1. We discuss the impacts of the uncertainty of the different datasets on the model calculation results.

3. Methods

3.1. Downscaling Method

(1)

Groundwater balance

In the study area, precipitation infiltration accounted for 86.3% of the total groundwater recharge based on the groundwater budget analysis [52,59]. Groundwater discharge is primarily caused by groundwater withdrawal and evapotranspiration (ET). Groundwater withdrawal has become a major method for groundwater discharge, supplying water for agricultural, municipal, industrial, and ecological uses [59,60], with more than 70% of the exploited groundwater being used for agricultural irrigation [61]. In the study area, most of the irrigation water is lost through ET, and the irrigation return flow is considered internalized water redistribution, which is implicitly accounted for in Equation (1):

P = ET + R + ∆TWS

(1)

where P is the precipitation, ET is the actual evapotranspiration, R is the net runoff, and ΔTWS is the total water storage change, including soil and groundwater storage changes.

Considering the scarcity of surface water in the arid area of the BTH region, the contribution of the surface water storage capacity, including rivers, lakes, and reservoirs, is neglected [7]. The change in the long-term average soil storage is relatively small in the NCP, and the water storage change is primarily driven by the groundwater storage change [62]. Therefore, the changes in groundwater storage are mainly subject to changes in precipitation and ET:

∆GWS = P − ET + r

(2)

where ∆GWS is groundwater storage change and r is the error due to the simplification.

(2)

Groundwater storage model equations

The aquifer system of the study area was divided into rectangular cells, the spatial resolution of which can be flexible: either 3°, the same size as the output from the GRACE data, or 0.05° or smaller. Based on Darcy’s law and the water balance principle, the groundwater storage changes in a certain cell should be equal to the lateral flux changes from its neighboring cells and the vertical flux changes from precipitation infiltrating and actual evapotranspiration, which leads to

$$\sum _i{T}_{ei}\frac{h_{ei}^{n+1}-h_e^{n+1}}{d_{ei1}+d_{ei2}}{L}_{ei}+A_e\sum _e{Q}^n=S_{ye}\frac{h_e^{n+1}-h_e^n}{\Delta t}{A}_e i=1,2,3,4$$

(3)

$$\sum _e{Q}^n={\alpha }_e{Q}_p^n-{\beta }_e{Q}_w^n+Q_c^n$$

(4)

where i is the number of cells around cell e; n is the nth time step; e_i represents the ith cell around cell e; $T_{ei}$ is the transmissivity for ith cell around cell e [L²T⁻¹]; $h_e^n$ and $h_{ei}^n$ are the groundwater level for cell e and its neighboring cell $ei$ at the nth time step, respectively [L]; $d_{ei1}$and $d_{ei2}$are the distance from the center of the cell and that of the cell neighboring from the adjacent boundary, respectively [L]; $L_{ei}$ is the lateral width of the aquifer [L]; $A_e$ is the area of cell e [L²]; ${\alpha }_e$ is the infiltration rate of precipitation (dimensionless); $Q_p^n$ is the monthly precipitation [LT⁻¹]; ${\beta }_e$ is the evaporation rate of the groundwater (dimensionless); $Q_w^n$ is the monthly evapotranspiration [LT⁻¹]; $Q_c^n$ is additional monthly groundwater recharge and discharge from such activities as local pumping [LT⁻¹]; $S_{ye}$ is the comprehensive specific yield (dimensionless); and $\Delta t$ is a time step chosen to be 1 month, the same as the time interval of the GRACE data [T].

The GRACE-derived groundwater storage anomaly is defined as the difference between the indicated groundwater storage and the average groundwater storage from January 2004 to December 2009 [L], which is denoted as (GW_g) in the equation. It is assumed that GWSA is accurate and error-free. The relationship between GWL and GW_g at a certain cell e can be expressed as:

$$G{W}_{ge}^n=S_y{}_e(h_{{}_e}^n-h_{avge})$$

(5)

where $h_{avge}$ is the average GWL for cell e from January 2004 to December 2009 [L].

Hence, Equation (3) can be rewritten as the following equations:

$$G{W}_{ge}^{n+1}(-\sum _i\frac{T_{ei}{L}_{ei}}{S_{ye}(d_{ei1}+d_{ei2})}-\frac{A_e}{\Delta t})=-\sum _i\frac{G{W}_{gei}^{n+1}{T}_{ei}{L}_{ei}}{S_{yei}(d_{ei1}+d_{ei2})}-\sum _i{T}_{ei}\frac{h_{avgei}-h_{avge}}{d_{ei1}+d_{ei2}}{L}_{ei}-A_e\sum _e{Q}^n-\frac{G{W}_{ge}^n}{\Delta t}{A}_e$$

(6)

It should be noted that observation wells are very limited and unevenly distributed in most conditions and that the initial hydraulic gradient in groundwater systems is difficult to achieve. Therefore, it is assumed that initial hydraulic gradient can be estimated by multiplying the surface slope by an unknown coefficient, which is denoted as hydraulic gradient coefficient.

$$T_{ei}\frac{h_{avgei}-h_{avge}}{d_{ei1}+d_{ei2}}{L}_{ei}\approx T_{ei}\frac{Z_{avgei}-Z_{avge}}{d_{ei1}+d_{ei2}}{C}_h{L}_{ei}$$

(7)

where $Z_{avge}$ and $Z_{avgei}$ are the average surface elevations at cell e and the neighboring cell e_i, respectively, and $C_h$ is the hydraulic gradient coefficient.

For each grid cell, we obtained a series of equations that may be resolved by considering the boundaries and initial conditions. Thus, the GW_g in each cell can be obtained from the equation. The principle of the model has also been applied in the Shiyang River basin in northwest China [63].

(3)

Downscaling

The downscaling can be divided into two steps. The first step is the groundwater storage model construction, calibration, and verification. The study area was divided into 1° × 1° grids in the horizontal direction. The number of simulation layers was one, and it was divided into 81 grids in total (Figure 3). The GWSA for January 2003 was set as the initial condition. The cells along all the boundaries (G50–G81) and cells G28, G34, G35, and G42 (which are close to the sea) are represented by the Dirichlet boundary with the GRACE-derived GWSA. The precipitation and actual evapotranspiration data were used to develop the model. According to regional groundwater studies in China [64], four types of aquifers are present in the study area (Figure 3). The parameters of the model included the comprehensive specific yield ($S_y$), transmissivity ($T$), precipitation infiltration coefficient ($\alpha$), evapotranspiration coefficient ($\beta$), and hydraulic gradient coefficient ($C_h$), which were determined according to empirical values from the different aquifers. The Morris one-at-a-time (MOAT) method was used for parameter calibration owing to its simplicity and efficiency [65,66]. The period from January 2003 to June 2017 was set as the calibration period, and June 2018 to December 2020 was set as the verification period (without considering the gap between GRACE and GRACE-FO). The time step was selected to be one month, which is consistent with GRACE data. The fitting objects of the model were the GRACE-derived time-series GWSA in 45 grid cells, except for the Dirichlet boundary grid cells. The groundwater storage model code (NGFLOW-GRACE) was used to calculate the GWSA for the study area.

The second step was to run the model with a refined mesh (with a spatial resolution of 0.05°) to obtain a higher-resolution GWSA. During the two processes, the distribution of the hydrogeological parameters remained the same and was only estimated in the first step, further spreading to 0.05° space. The model parameters for the final calibration are discussed in Section 4.3. The precipitation and evapotranspiration data with a resolution of 0.05° were used to drive the model and obtain the GWSA at small scales.

3.2. Evaluation Methods

The GRACE-derived GWSA was marked as $G{W}_{g\_obs}$, and the GWSA simulated from the NGFLOW-GRACE model was denoted as $G{W}_{g\_sim}$. The model performance was measured using total root mean square error (Φ), as shown below:

$$\Phi =\sqrt{\frac{1}{N}{\sum }_{m=1}^N{\sum }_{n=1}^{Ndt}{\left(G{W}_{mg\_sim}^n-G{W}_{mg\_obs}^n\right)}^2}$$

(8)

where N is the number of cells with an unknown GWSA in the model area. The time step is Ndt, and m and n are the continuous counts of grid cells and time steps, respectively.

A shuffled complex evolution method was used for parameter estimation [65]. The parameter optimization strategy is to minimize Φ.

The root mean squared error (RMSE) and Nash–Sutcliffe efficiency coefficient (NSE) were used to assess the performance of the model data products. The RMSE was used to further measure the deviation between the simulated and GRACE-derived and ground-based observations in each grid cell. The NSE was used to evaluate the degree of fit of the regression models. An RMSE close to 0 and NSE close to 1.0 indicate that the model performance is reliable.

$$RMSE=\sqrt{\frac{1}{n}{\sum }_1^n{\left(G_i-P_i\right)}^2}$$

(9)

$$NSE=1-\frac{{\sum }_{i=1}^n{\left(G_i-P_i\right)}^2}{{\sum }_{i=1}^n{\left(G_i-\overline{G_{𝜄}}\right)}^2}$$

(10)

where n is the total amount of data, $G_i$ is the GRACE-derived results or ground-based observations, $P_i$ is the model-simulated results, and $\overline{G_{𝜄}}$is the average value.

4. Results

4.1. Model Evaluation

The calibration period of the model was from January 2003 to June 2017. The TRMM (precipitation) and GLEAM (AET) data were first used to drive the model, in which TRMM data have been proven to match the in situ observations in the study area well [7]. Figure 4 illustrates that the GWSA calculated by the model and the GWSA by GRACE-derived data are consistent during the calibration period. Figure 5 shows the model evaluation metrics RMSEs and NSEs for the 45 grid cells, where the RMSEs change from 0.81 cm EWH to 5 cm EWH with an average value of 2.61 cm EWH and the NSEs vary from 0.41 to 0.96 with an average value of 0.84. During the validation period from June 2018 to December 2020, the trend of the simulated values and GRACE-derived values remain similar, with a minimum RMSE of 1 cm EWH and a maximum NSE of 0.83. Overall, the simulated results agree well with the GRACE-derived results. The simulated GWSA matched the GRACE-derived results the best for the northwestern mountainous area, where the changes in GWSA were less affected by human activities.

4.2. Downscaled GWS Changes

After estimating the hydrogeological parameters, the model was refined from a spatial resolution of 1° (Figure S1a) to 0.05° (Figure S1b). Therefore, a grid cell at a spatial scale of 1° was refined to 400 sub-grid cells at a spatial scale of 0.05°, and the spatial parameters in Figure 3 were averaged from a resolution of 1° to 0.05°. For the model with a spatial scale of 0.05°, the PENG precipitation and MOD16 AET data were used to drive the model. The initial conditions and boundary conditions were the same as those used in the calibrated model. The simulation period for the downscaled models was from January 2003 to December 2020 without considering the gap between the GRACE and GRACE-FO data.

To better demonstrate the characteristics of the downscaled results, the GWSA of a 1° grid cell was compared with that with a spatial resolution of 0.05° for each of the 45 grid cells (Figure S2). The downscaled GWSA maintained the same pattern as that of the same 1° grid. The correlation coefficient of the 0.05° average GWSA within each grid and the GRACE-derived GWSA in the same 1° grid were above 0.9. Figure 6 contrasts the GRACE-derived GWSA variations with the downscaled GWSA variations on a spatial scale for 2019. Three transects (P1–P2, P3–P4, and P5–P6 across 4°, 4°, and 1°, respectively) demonstrate the high-resolution GWSA data along the transect capture the spatial variations of the GWSA, which appear to be uniform in the coarse-resolution GRACE data. High-resolution GWSA data are particularly useful for small regions in which the downscaled GWS data capture the distribution variations. For example, the downscaled data reflect a strong spatial heterogeneity in the GWSA variations (−22.98 to −39.05 cm EWH) in transect P5–P6 with a substantial decline observed, whereas the GRACE-derived GWSA is represented by homogenized pixels. Therefore, the downscaling results not only exhibit a time series variation trend consistent with that of the GRACE-derived values but also exhibit finer spatial heterogeneity.

4.3. Validation of Downscaling Results

The time series of the GRACE-derived and downscaled GWSA are compared with the in situ GWL data in Figure 7. The result showed that the GWSA agrees well with the in situ GWL observations, that is, the downscaled results were also capable of capturing the significant groundwater change information in most cases. The GWSA had a significant downward trend from 2005 to 2014, which can be attributed to excessive demand for groundwater required for agriculture during the 2003–2009 precipitation deficiency period and average annual groundwater withdrawal of 2.5 × 10¹⁰ m³, resulting in unsustainable recovery of groundwater storage [7,11].

Figure 8 quantifies the comparison between the GRACE-derived, downscaled GWSA and in situ observed GWL in terms of correlation coefficients (CC) and RMSE. For the computation of RMSE and CC, two time series were normalized to [−1,1] because their units were inconsistent. Downscaling the GRACE data led to an increase in the CC and a decrease in the RMSE, indicating the downscaling method improved GWSA computation. On average, higher correlations are observed after downscaling, which are improved from 0.78 to 0.85. The best cases occur for the H52 well, with the CC increasing from 0.74 to 0.91 and the RMSE decreasing from 0.34 to 0.21 in the B4 well. Overall, the results indicate that the downscaled process increases the correlation coefficients between the GRACE-derived GWSA and in situ GWL measurements by 0.06 (on average) and reduces the RMSE by 8.85% (on average). The least agreement between the downscaled results and GWL was found with the shallow monitoring wells. The key reasons for this could be due to the exploitation of the deep groundwater.

5. Discussion

5.1. Spatiotemporal Changes of Downscaled GWSA in Subregions

The linear least squares regression method was used to calculate the trends of the variables in the current study. Figure 9 presents the temporal fluctuations of the monthly GWSA in the BTH region and three subregions (PP, ECP, and WM) from January 2003 to December 2020. The GWS declines over the BTH with different decline rates during the four periods, i.e., January 2003–December 2008 (−0.06 cm/month), January 2009–December 2014 (−0.097 cm/month), January 2015–June 2017 (−0.216 cm/month), June 2018–December 2020 (−0.143 cm/month). The time series demonstrate that the three subregions exhibit similar variations; GWSA has a decreasing trend, ranging from −0.053 to −0.071 cm/month in January 2003– December 2008, followed by a larger fluctuation and decrease rate in the PP from January 2009 to December 2014. However, GWS has a continuously accelerated decline with a large difference in the subregions from January 2015 to June 2017, and the rate of decrease is more than two times those of the previous period, which are −0.313 cm/month, −0.183 cm/month, −0.239 cm/month, and −0.216 cm/month in the PP, WM, ECP, and BTH, respectively; similar trends have been reported in the literature [67,68]. The western mountainous areas showed the slowest decline in the different periods. Since 2018, the rate of decline has slowed in the BTH and the three subregions.

Figure 10 shows the spatial patterns of the linear trends of the 0.05° GWSA with the same legend range from −0.562 to −0.084 cm/month. For the spatial distribution, the GWSA generally declined from 2003 to 2020, while areas of groundwater storage decline are more visible in the PP. The downward trend gradually increased from northeast to southwest, and the maximum decline rate appeared in Handan City, Hebei Province. From January 2015 to June 2017, the decline rate of the GWSA in the PP was the most serious and was mainly distributed in Handan, central Xingtai, eastern Shijiazhuang, and western Cangzhou. In addition, the decline in the GWSA in the Zhangjiakou area intensified during this period. From June 2018 to December 2020, especially in Beijing, the GWS decline slowed. The coastal areas south of Tianjin and west of Cangzhou experienced stable GWS changes over the years. In the BTH region, the decline rate of the GWS gradually decreased owing to the groundwater overexploitation control policy in recent years.

5.2. Comparison of GWS with GWL Observations in Administrative Regions

To demonstrate the relationship between the GWL changes and GWS changes in the study area, the downscaled GWS changes were compared using the GWL changes from 843 observation wells in the subregions, including Cangzhou, Xingtai, Handan, Tangshan in Hebei Province, Beijing, and Tianjin. Because of the different time series of the collected monitoring well data, the GWL and GWS data for Cangzhou, Xingtai, Handan, and Tangshan from December 2018 to December 2019 were selected for comparison. The monitoring wells were divided into shallow and deep wells. Beijing and Tianjin monitoring wells from December 2014 to December 2016 were shallow wells. The GWS was not represented by layers (because the GWS derived from GRACE was difficult to divide vertically). It should be noted that owing to the true value of specific yield is hard to obtain, we only compared groundwater level data from observation wells with groundwater storage changes in the pixel including the wells. The GWS changes legend for Beijing and Tianjin were different from those in the other three subregions.

Figure 11 shows that the change trend of the deep monitoring wells and that of the GWS changes in the spatial location correspond well; in other words, at locations where the GWL drops more, the GWS decreases more, such as at the locations in the red dotted circle. However, there are also inconsistencies, such as the locations in the blue dotted circle; the GWL drops are the largest, but the GWS decline is relatively small. This could be attributed to the local groundwater extraction, which is not considered in the model due to the lack of detailed information.

Figure 12 shows that the GWS and GWL changes also exhibited consistent trends in the monthly variations. As shown in Figure 12a–d,f, from March to June, the GWS gradually decreased and the GWL also gradually decreased, mainly due to groundwater exploitation for agricultural irrigation. With the increase in precipitation in July, the GWS and GWL showed an increasing trend. In August, the GWL demonstrated an upward trend, but the GWS continued to decrease, which is difficult to explain. From October to December, when groundwater exploitation was slow, the GWS gradually recovered and the GWL also slowly increased. However, Figure 12e shows that the variation trends of GWS and GWL are inconsistent, which indicates that in the absence of deep-monitoring GWL verifications, the variation trend of the GWS cannot adequately explain the variation in the shallow GWL in Beijing.

Table 2. Statistics of GWL observations and GWS changes in the subregions.

Time	December 2018–December 2019				December 2014–December 2016
Sub-Regions	Cang Zhou	Xing Tai	Han Dan	Tang Shan	Tian Jin	Bei Jing
Number of shallow wells (up)	16	21	4	7	17	69
Number of shallow wells (down)	118	90	48	23	26	67
Number of deep wells (up)	31	13	1	----	----	----
Number of deep wells (down)	145	117	30	----	----	----
GWL changes in shallow aquifers (m)	−0.71	−1.15	−2.36	−0.42	−0.98	0.29
GWL changes in deep aquifers (m)	−2.30	−3.63	−3.76	----	----	----
GWS changes (cm EWH)	−1.37	−2.55	−3.65	−1.21	−5.06	−6.14

summarizes the GWL changes in the monitoring wells and GWS changes in each subregion during the different time periods. Between December 2018 and December 2019, the GWL in the shallow aquifers of Tangshan, Cangzhou, Xingtai, and Handan dropped by 0.42 m, 0.71 m, 1.15 m, and 2.36 m, respectively, and the average GWL in the deep aquifers dropped by 2.30 m, 3.63 m, and 3.76 m (without the Tangshan), respectively. The decrease in the GWL in the deep aquifers was generally greater than that in the shallow aquifers, indicating that the main groundwater consumption was from the deep confined aquifer. The GWS decreased on average by 1.21 cm EWH, 1.37 cm EWH, 2.55 cm EWH, and 3.65 cm EWH, respectively. The good correspondence between the GWL and GWS changes indicated that groundwater consumption was the most serious in Handan.

From December 2014 to December 2016, the GWS in Beijing and Tianjin decreased by 6.14 cm EWH and 5.06 cm EWH, respectively, and the GWS in southwest Beijing and central Tianjin declined the most. From the perspective of GWL changes, the average GWL in Tianjin dropped by 0.98 m, while the average GWL in Beijing rose by 0.29 m, which is closer to the 0.43 m increase from 2014 to 2016 listed in the Beijing Water Statistical Yearbook in 2016 [69]. The GWL and GWS changes in Beijing exhibited completely different trends, which were similar to the results of GRACE-derived GWS changes and GWS changes estimated using GWL by Zhang et al. [50]. There may be two reasons for this result. First, the uncertainty of the GRACE data because the subregion GRACE data are contaminated by signal leakage and amplitude-damping effects [11,22,50]. Second, the water from the South-to-North Water Diversion began to enter Beijing in December 2014, replacing the exploitation of some groundwater sources, and the surplus water of the South-to-North Water Diversion was used for groundwater replenishment. The GWL in Beijing has gradually recovered since 2015; however, the effects related to human factors were not considered in the model.

5.3. Evaluation of Hydrogeological Parameters

Figure S3 presents the spatial distribution of the relevant parameters in the downscaled model. The comprehensive specific yield ranged from 8.0 × 10⁻⁴ to 1.0 × 10⁻³ in the piedmont plain gradually decreased to 4.0 × 10⁻⁵ to 1.0 × 10⁻⁴ in the coastal areas, and the intermountain basins had the smallest comprehensive specific yield of less than 4.0 × 10⁻⁵. This paper compares the parameters results with previous studies. Previous research has shown that the specific yield of shallow aquifers in the plain area of the study area is between 0.025 and 0.29, the storage coefficient in the deep aquifers decreases from 1.0 × 10⁻³ to 8.0 × 10⁻³ in alluvial plains to 4.0 × 10⁻⁴ to 5.0 × 10⁻⁴ in coastal plains, and the average storage coefficient of deep aquifers is 0.00125, approximately 2% of that of the shallow aquifer area-weighted specific yield (0.075) [22,47,50,59]. Shao et al. [59] determined that the specific yield of the aquifer exhibits a downward trend from the piedmont alluvial plain to both sides of the river valley, that the specific yield ranges from 0.03 to 0.23, and that the order of magnitude of the storage coefficient in the confined layer is 10⁻⁴–10⁻⁶. Li et al. [52] calculated specific yield ranging from 0.025 to 0.2 and storage coefficients in the confined layer ranging from 2.7 × 10⁻⁵ to 2.3 × 10⁻⁴, which were provided in the numerical simulation. The model simulated the comprehensive changes of phreatic and confined aquifer. Thus, we redefined traditional specific yield as comprehensive specific yield in this study, as mentioned above. We found that our optimized storage coefficient has the same order of magnitude as the elastic storage coefficient of confined water. The groundwater exploitation in the study area is dominated by deep confined water. Therefore, we think the optimized comprehensive water supply is reasonable.

The precipitation infiltration coefficients in the study area ranged from 0.02 to 0.30. The precipitation infiltration coefficients were greater than 0.2 for the piedmont plains, between 0.05 and 0.2 for the central and eastern plains, and less than 0.03 for the mountainous and coastal plains. This result is similar to the parameters of the plain area optimized by Li et al. [52], which was spatially distributed within the range of 0.17 to 0.25. The evapotranspiration coefficient ranged from 0.03 to 0.66. The evapotranspiration coefficients were greater than 0.50 in the piedmont plains—less than 0.1 in most areas of the intermountain basin, 0.1–0.25 in the coastal areas, and 0.25–0.5 in most of the central and eastern plains, which is consistent with the parameters in the Handbook of Hydrogeology [70].

5.4. Uncertainty Analysis

The accuracy of a model is restricted by uncertainties in climate forcing (particularly precipitation and AET) and model parameters as well as deficiencies in the model structure. The MOAT method with 1000 samplings for each parameter was used to determine the sensitivity of the model to the aquifer parameters. As shown in Figure S4, the hydraulic gradient coefficient and comprehensive specific yield for the porous aquifer had the highest degrees of sensitivity. The parameter with the lowest degree of sensitivity was transmissivity. Given that precipitation and evapotranspiration are the main forcing data of the model, 12 scenarios with different precipitation and AET data inputs were set to examine the influences of the model input data. Figure S5 illustrates that the GWSA trends from the twelve scenarios exhibit the same patterns over the entire area, and the maximum absolute difference in the GWS changes is less than 1 cm EWH, indicating that although there are some differences between precipitation and evapotranspiration from different data sources, the different forcing data may not have substantially influenced the groundwater storage changes. Therefore, the average groundwater storage values for the different forcing data were not further implemented.

5.5. Advantages and Limitations

The groundwater storage model was established using the free GRACE data applied to improve the GWSA resolution from 1° to 0.05° and assess the groundwater storage changes at small scales in this study. This study provides a new method for the downscale study of GRACE data, which is based on the hydrological process of groundwater. This is an example of the combination of remote sensing, in situ monitoring, and hydrological models. Remote sensing data provide driving data and constraints for the model, and in situ monitoring provides an important data source for model verification. This study can provide a new idea for the combination of remote sensing data and numerical models and form an important dataset of groundwater storage changes for water resource management. However, the present study has some limitations. First, the established model considered the GRACE-derived GWSA as the fitting target; therefore, the accuracy of the GRACE-derived GWSA data also affects the accuracy of the downscaling results [71]. In this model, it is assumed that the GRACE-derived GWSA data are error free. Second, after generalizing the water balance in the study area, this study only considered groundwater storage changes driven by natural factors, such as precipitation and actual evapotranspiration. Surface water changes, such as reservoir scheduling, water transfer projects [53], and local groundwater extraction, were not fully considered in the model. Although the research results are relatively reasonable, to provide more comprehensive scientific guidance for policy making, the impact of human activities on groundwater storage changes must be considered in the model [72]. Third, comprehensive specific yield and hydraulic gradient coefficient were the most sensible variables. Whether those sensitive parameters will change with time and downscaled spatial is not discussed in this study. Furthermore, the observational data used for the evaluation were limited, and those will be the focus of future research.

6. Conclusions

Because the use of GRACE products for smaller areas is limited due to its coarse spatial resolution, this study utilizes a dynamic downscaling method to improve the GWSA resolution from 1° to 0.05° by constructing a groundwater storage numerical model with GWSA as the key variable. First, we constructed groundwater storage model and performed model calibration and validation and verified the downscaling results. Then, the downscaling results were used to evaluate the temporal and spatial changes in groundwater storage at both the regional and local scales and were compared to the groundwater level changes in different administrative regions. In addition, the model parameters and uncertainty analysis were discussed. This study demonstrates the potential of downscaling GWSA data for local-scale applications. The following conclusions were drawn from this study:

(1) The study area is divided into 1° × 1° grids with one layer in the horizontal direction, and the initial GWSA from January 2003, Dirichlet boundary, and initial parameters are given. After parameter optimization, the simulated GWSA at the spatial resolution of 1° had an average RMSE of 2.61 cm EWH and average NSE of 0.84 for the calibration period and the minimum RMSE of 1.0 cm EWH and maximum NSE of 0.83 for the validation period. The calibrated parameter distributions are consistent with those from the field investigation.

(2) The downscaled results exhibited a time series variation trend that was consistent with that of the GRACE-derived results, of which the correlation coefficient between the average GWSA within 0.05° grid and the GRACE-derived GWSA in the same 1° grid was above 0.9. The results also showed finer spatial heterogeneity. Compared to the GRACE-derived GWSA, the downscaled process increased the correlation coefficients between the GWSA and in situ measurements of the GWL by 0.06 (on average) and reduced the RMSE by 8.85% (on average).

(3) Based on the downscaling results, the variation pattern of groundwater storage demonstrated that the downward trend gradually increased from northeast to southwest in the study area, while groundwater storage in the piedmont plains declined significantly by 0.313 cm/month from January 2015 to June 2017, and the western mountainous areas declined the slowest during the different periods. However, the decline rate of the GWS has gradually decreased owing to the groundwater overexploitation control policy in recent years.

(4) The comparison between the GWL and GWS changes for different administrative regions showed good correspondence in both the spatial and monthly variations and indicated that groundwater consumption was the most serious in Handan, Hebei Province, where the GWL declined by 2.36 m and 3.76 m in the shallow and deep aquifers, respectively, and the GWS decreased by 3.65 cm EWH from December 2018 to December 2019, followed by Xingtai, Hebei Province.

(5) The MOAT method with 1000 samplings for each parameter in the model was used to carry out a sensitivity analysis, which indicated that the hydraulic gradient coefficient and comprehensive specific yield have the highest sensitivity and that transmissivity has the lowest sensitivity. In addition, different forcing data had no substantial influence on the groundwater storage changes.

The proposed method takes GRACE-derived groundwater storage data as the fitting variable of the model, which can be applied in remote areas where observation well data are lacking. Meanwhile, a dataset of groundwater storage changes with a higher resolution can be obtained and may be public for the further analysis of groundwater-related problems. For future work, the downscaling of groundwater storage changes at refined time scales will be conducted, as will the improvement of model reliability at refined spatial scales. Then, the relationship between downscaled groundwater storage changes and groundwater discharges under climate changes will further be explored to estimate the true groundwater withdrawals for effective groundwater management.