[go: up one dir, main page]
Next Article in Journal
Small Object Detection Based on Deep Learning for Remote Sensing: A Comprehensive Review
Previous Article in Journal
Near Real-Time Flood Mapping with Weakly Supervised Machine Learning
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 15
Issue 13
10.3390/rs15133264
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Reconstructing Groundwater Storage Changes in the North China Plain Using a Numerical Model and GRACE Data

by
Junchao Zhang
1,2,
Litang Hu
1,2,*,
Jianchong Sun
1,2 and
Dao Wang
1,2
1
College of Water Sciences, Beijing Normal University, Beijing 100875, China
2
Engineering Research Center of Groundwater Pollution Control and Remediation of Ministry of Education, Beijing Normal University, Beijing 100875, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(13), 3264; https://doi.org/10.3390/rs15133264
Submission received: 27 May 2023 / Revised: 22 June 2023 / Accepted: 24 June 2023 / Published: 25 June 2023
Browse Figures
Figure 1
<p>Study area of the NCP.</p> "> Figure 2
<p>Comparison of different precipitation and actual evapotranspiration data sources in the NCP from 2003 to 2019: (<b>a</b>) monthly precipitation; (<b>b</b>) actual evapotranspiration.</p> "> Figure 3
<p>Distribution of parameter zones and Dirichlet boundaries (<b>a</b>) and geological cross section of the NCP (<b>b</b>) in the study area.</p> "> Figure 4
<p>Simulated GWSA from January 2003 to December 2019.</p> "> Figure 5
<p>Comparison of groundwater storage changes before and after downscaling in December 2019. (<b>a</b>,<b>c</b>) The 1° and 0.05° mesh grid blocks, respectively; (<b>b</b>,<b>d</b>) the spatial distribution of GWSA before and after downscaling.</p> "> Figure 6
<p>Time-series trend of downscaling results from 2003 to 2020. The red bands represent the changes in GWSA within all 0.05° grids under the 1° grid, and the blue lines represent the GWSA retrieved by GRACE and GLDAS in the 1° grid.</p> "> Figure 7
<p>Comparison of time-series groundwater storage changes from the model and field observations.</p> "> Figure 8
<p>Comparison of TWSA of (<b>a</b>) Li data, (<b>b</b>) CLSM data, and (<b>c</b>) Mo data before and after standardization with the GRACE-retrieved data.</p> "> Figure 9
<p>Change in GWSA from 2003 to 2020 (the GWSA data from July 2017 to May 2018 were calculated using Li data, Mo data, and GLDAS CLSM model data).</p> "> Figure 10
<p>Change in the GWSA in the NCP and administrative regions (<b>a</b>) and change in GWSA and precipitation in the NCP (<b>b</b>) from 2003 to 2020.</p> "> Figure 11
<p>Spatial distribution map of the GWSA change rate in the NCP in four subperiods: (<b>a</b>) 2003–2008, (<b>b</b>) 2009–2014, (<b>c</b>) 2015–2017, and (<b>d</b>) 2018–2020.</p> "> Figure 12
<p>Spatial distribution map of the GGDI change rate: (<b>a</b>) 2003.1 to 2008.12, (<b>b</b>) 2009.1 to 2014.12, (<b>c</b>) 2015.1 to 2017.12, and (<b>d</b>) 2018.1 to 2020.12.</p> ">
Review Reports Versions Notes

Abstract

:
Groundwater has been extensively exploited in the North China Plain (NCP) since the 1970s, leading to various environmental issues. Numerous studies have utilized Gravity Recovery and Climate Experiment (GRACE) satellite data to analyze changes in groundwater storage in the NCP and provide valuable insights. However, the low spatial resolution of GRACE data has posed challenges for its widespread application, and there have been limited studies focusing on refining groundwater storage changes in the NCP. In addition, the lack of data on the gap period between GRACE and GRACE-FO hinders in-depth research on regional groundwater storage anomalies (GWSA). This paper applied a groundwater storage model called NGFLOW-GRACE to construct a groundwater storage change model in the NCP at spatial resolutions of both 1° and 0.05°. The groundwater storage change model was calibrated and driven using gratis data, with hydrogeological parameter values estimated using the shuffled complex evolution algorithm (SCE-UA). The model exhibited favorable performance, with correlation coefficients greater than 0.85 during the calibration period and 55% of coefficients greater than 0.50 during the validation period. Interestingly, the results indicate that different combinations of remote sensing data do not significantly impact the outcomes, while the hydraulic gradient coefficient demonstrates the highest sensitivity. Appropriate reconstructed data were selected within the empty window period, and by downscaling the model to a resolution of 0.05°, a complete cycle (January 2003 to December 2020) of GWSA was derived. Through comprehensive comparisons with previous research findings on both temporal and spatial scales, it can be concluded that the downscaled groundwater storage changes obtained from the established model demonstrated high reliability.

1. Introduction

As the main source of fresh water in many regions of the world, groundwater supplies approximately 50% of the demand for drinking water, 40% of the demand for self-sufficient industries, and 20% of the demand for irrigation water [1]. Furthermore, groundwater can not only store a large amount of precipitation but can also effectively maintain the flow of rivers in dry periods [2]. However, globally, the increases in extreme climate and the frequency of human activities such as groundwater overexploitation have led to aquifer depletion [3]. For example, in India, groundwater exploitation has been increasing at an alarming rate, leading to severe groundwater depletion in Punjab, Haryana, parts of Rajasthan, and Uttar Pradesh [4]. In the North China Plain (NCP), excessive groundwater exploitation has led to groundwater depletion and pollution [5]. In Mexico, population growth and industrial development in some regions have led to increased water extraction to the point where future water supplies are threatened [6]. The depletion of aquifers not only limits the sustainable development of the regional economy but also brings about a series of serious ecological and environmental problems, which are regarded as an acute threat to the sustainability of the water supply [7].
The excessive exploitation of groundwater in the NCP is of great concern. In 2013, the state determined the water receiving areas (WARs) of the South to North Water Transfer (SNWTP) as Beijing, Tianjin, Hebei, Shandong, and Henan in an effort to replace urban groundwater exploitation through the use of the SNWTP. In December 2014, the central line of the SNWTP, which starts from the Danjiangkou Reservoir and ends at the Summer Palace in Beijing, was determined to be open to water, with beneficiary areas including Henan, Hebei, Beijing, and Tianjin. In 2019, the Action Plan for Comprehensive Treatment of Groundwater Overexploitation in North China was issued and implemented. The plan systematically promotes the treatment of groundwater overdraft in North China, aiming to gradually realize the balance of groundwater extraction and replenishment. In addition, the first groundwater management administrative regulation, Groundwater Management Regulations, was promulgated and implemented in December 2021. Significant results have been achieved regarding the treatment of groundwater overdraft in the NCP, but further study is still required to improve the management of groundwater and strengthen the monitoring of groundwater resources.
Point-based measurement, satellite-based monitoring, and regional groundwater estimation through numerical modeling are the most original methods and tools for detecting groundwater level fluctuations [2]. Point-based measurements take time and manpower and require a great number of monitoring points to ensure their representativeness. The accurate quantification of large-scale regional groundwater storage changes is often limited by various factors, and relying solely on field groundwater measurements can pose challenges, particularly when trying to capture long-term changes on a regional scale. Hydrological models require more comprehensive meteorological and hydrological information and tend to ignore the complexity of key hydrological processes, resulting in inaccurate detection results. Consequently, the shortcomings of traditional groundwater monitoring methods have been unable to meet the current needs for the management of groundwater resources. The Gravity Recovery and Climate Experiment (GRACE) twin satellites made it possible to explore groundwater storage (GWS) at regional scales. Additionally, GRACE Follow-On (GRACE-FO), which can sense the changes in the Earth’s gravity field at a smaller scale and realize accurate observations of the Earth’s gravity, was launched in 2018 to replace the GRACE satellite [8]. Although the spatial and temporal resolution (no more than 160,000 km2 weekly or monthly) of GRACE is relatively low compared with other Earth observation satellites, the prime advantage of GRACE is that it can detect water stored at all levels, including groundwater [9]. If other water storage components, such as soil moisture, surface water, and ice water, can be accurately estimated from the model or observation and removed, the groundwater storage changes could be estimated from the terrestrial water storage changes derived from GRACE [10]. In the past decade, research on groundwater storage estimation related to GRACE has achieved fruitful results and has been applied worldwide, e.g., India [11], the NCP [5], and the Central Valley of the United States [12]. Huang et al. [13] used GRACE data to discuss the shallow and deep groundwater storage changes in the NCP, and the estimation of shallow groundwater storage changes had the highest correlation with the field-measured data, with the correlation coefficient (R2) reaching 0.91; additionally, the R2 of the estimation of the deep groundwater storage changes and the field-measured data reached 0.75. Many previous studies have confirmed that the groundwater storage in the NCP has been declining [5,14,15]. However, the low resolution of GRACE data, which cannot be fundamentally solved by data postprocessing, is still an intractable difficulty for the application of GWS change worldwide [15]. Consequently, drawing upon downscaling methods to improve data spatial resolution is essential to ultimately contribute to regional groundwater management.
Downscaling, effectively, is an approach used to obtain high-resolution data from low-resolution data. In general, the statistical downscaling method and dynamic downscaling method are the two most commonly used downscaling methods. The statistical downscaling method establishes statistical relationships between large-scale climatic characteristics and regional climatic elements using years of observational data [16]. For example, Yin et al. [17] used evapotranspiration data to establish a statistical relationship with the groundwater storage anomalies obtained by GRACE RL05 inversion, which reduced groundwater storage results retrieved by GRACE from 110 km to 2 km in the NCP. Rahaman et al. [18] successfully reduced the groundwater storage change result retrieved by GRACE RL05 from 1° × 1° to 0.25° × 0.25° using the random forests method. Seyoum et al. [19] used the boosting regression tree method to establish the regression statistical relationship with the terrestrial water storage anomaly (TWSA) from GRACE, and obtained high-resolution groundwater level anomalies.
Although the statistical downscaling method has the advantages of less calculation and relatively uncomplicated model construction, it lacks a strong physical basis and is not applicable for areas where the correlation between large-scale features and regional elements is weak. As opposed to the statistical downscaling method, dynamic downscaling has a solid foundation in mathematical physics and is based on regional climate models. Global models provide initial boundary conditions for regional climate models and obtain high-resolution weather and climate information through numerical integration of high-resolution regional climate models [20]. The data assimilation processing of the numerical model based on the terrestrial water storage changes estimated by GRACE satellite data can effectively reduce the uncertainty of the regional water storage changes estimated by the numerical model and improve the reliability of the simulation. For example, Houborg et al. [21] improved the GRACE assimilation model proposed by Zaitchik [22] and used GRACE RL04 as the research data to apply the assimilation model to drought monitoring in North America. The research results highlighted the application potential of the GRACE assimilation model in drought monitoring. Previous studies have shown the applicability of this model in different regions [23,24]. Therefore, the focus of this article was to use reconstructed data to fill the gap between GRACE and GRACE-FO data. By conducting a thorough comparative analysis of the reconstructed data from previous studies, our objective was to identify and employ more suitable window data for the downscaling process. This approach aimed to enhance the accuracy and reliability of our downscaling method and enabled a comprehensive spatiotemporal analysis of the downscaled results.
The objectives of this study were to (a) establish a groundwater storage model with a 1° spatial resolution using a dynamical downscaling model named the NGFLOW-GRACE model in the NCP from 2003 to 2020 and conduct calibration and validation of the model to obtain hydrogeological parameter estimations; (b) establish a 0.05° model based on a 1° model and validate the downscaled GWSA with in situ observation data; and (c) compare the reconstructed data from previous studies and select the most suitable window data to obtain a more precise and comprehensive comprehension of the spatial and temporal variations in groundwater storage in the NCP and five administrative regions from 2003 to 2020. The rest of this paper is organized as follows: Section 2 briefs the generality of the NCP. Section 3 describes the methodology of the model and data resources. Section 4 describes the establishment of the model in the NCP and the results after downscaling. Section 5 presents the results of the spatial distribution and time series after downscaling and a comparison of the downscaled results with previous results, and discusses the time series map and spatial distribution map of the GRACE groundwater storage drought index. Section 6 presents the conclusions.

2. Study Area

The NCP refers to the plain area east of the Taihang Mountains, south of the Yanshan Mountains, north of the Yellow River and east of the Bohai Sea, with geographic coordinates of 112.5°–119.5°E and 34.8°–40.4°N, that is, the plain area of the Haihe River Basin, covering an area of approximately 140,000 km2, and including 21 prefecture-level cities such as Beijing (BJ), Tianjin (TJ), Hebei Plain (HB), and North Henan (HN) and North Shandong (SD) areas north of the Yellow River (Figure 1). Based on major figures from the 2020 population census of China, the population in the NCP is approximately 110 million.
The surface elevation of the NCP does not exceed 100 m, and the terrain inclines toward Bohai Bay from the north, west and south directions. According to their genesis and morphological characteristics, they can be divided into piedmont alluvial plains, central and eastern alluvial plains, and alluvial and coastal plains. The NCP is formed by the sediment deposits of the Yellow River and the Haihe River, including the Yellow River, Luanhe River, Haihe River, the rivers in northern Shandong, and Tumajia River, as well as approximately 60 other rivers. The NCP has a semiarid and semi humid continental monsoon climate. Affected by the East Asian monsoon, the precipitation decreases from southeast to northwest. Approximately 70–80% of the precipitation occurs in the monsoon season from June to September. The average annual precipitation ranges from 500 mm to 600 mm, and the annual average temperature is between 8 °C and 15 °C. The annual average temperature and annual precipitation decrease with increasing latitude from south to north. The annual average potential evaporation is between 1100 and 2000 mm [25]. The NCP is likely to be arid because of its semi humid to semiarid climate, which highlights the need for the sustainable management of water resources in the region.
The Quaternary aquifer medium in the NCP is mainly composed of multiple types of superimposed aquifer group structures, with extremely complex and varied geometric forms, where the sand, gravel layers, and clay layers between them overlap and crisscross, forming the pore aquifer group [26]. As one of the largest aquifer systems in the world, the Quaternary aquifer in the NCP has tremendous groundwater resources and provides inestimable support for the country’s food production, agricultural employment, and industrial water delivery [27]. In light of the lithological distribution and hydraulic characteristics of the sediments in the NCP and in combination with the development and utilization of groundwater, the aquifer can be divided into four aquifer groups. The first aquifer group is an unconfined aquifer, which is the main water supply aquifer with a buried depth of 40–60 m, and the 2nd, 3rd, and 4th aquifer groups are confined aquifer groups with burial depths of 120–170 m, 250–350 m, and 350–550 m, respectively, and they are comprised of sand, gravel, clay, and silt [5].
The main type of recharge of shallow groundwater in the NCP is precipitation infiltration. However, in the past two decades, with the rapid increase in irrigation water and the development of urbanization, the evaporation of shallow groundwater has been significantly reduced, and groundwater exploitation has become the main groundwater discharge item. The deep groundwater is predominantly recharged by the hidden carbonate karst water in the piedmont edge and the lateral runoff and vertical infiltration recharge of the main alluvial fan in the piedmont. In Northern China, the groundwater supply accounts for nearly 30% of the total water supply [28]. The consumption of groundwater resources in the NCP is in an acute state. The use of groundwater storage in the NCP accounts for 50% to 75% of the total water consumption in the Huaihe River and Haihe River Basins [29]. Agricultural irrigation is the main cause of serious groundwater exploitation [30]. In the winter and spring dry seasons, the groundwater in the NCP maintains approximately 70% of irrigation and is mainly used for wheat production [31].

3. Methods and Data

3.1. Groundwater Storage Model

Based on the traditional groundwater mathematical model, the model (NGFLOW-GRACE) used in this study deduced the mathematical expression formula of groundwater storage changes according to the Darcy formula and water balance principle. In a certain cell, where the size of a cell is 1° or 0.05°, changes in groundwater storage are the total of the changes in water volume flowing into the cell laterally from its adjacent cells and the changes in vertical flux of precipitation infiltration, actual evapotranspiration, and other sources and sinks [23]:
∑ e T e i h e i n + 1 − h e n + 1 d e i 1 + d e i 2 L e i + A e ∑ e Q n = S y e h e n + 1 − h e n ∆ t A e   i = 1 , 2 , 3 , 4
∑ e Q n = α e Q p n − β e Q w n + Q c n
where T e i is the transmissivity for cell e [L2T−1]; h e n and h e i n are the groundwater level for cell e and its neighboring cell e i at the n t h time step, respectively [L]; d e i 1 and d e i 2 are the distance from the center of the cell and that of the cell neighboring the adjacent boundary, respectively [L]; L e i is the lateral width of the aquifer [L]; A e is the area of cell e [L2]; Q n is the vertical source and sink term [LT−1]; α e is the infiltration rate of precipitation [dimensionless]; Q p n is the monthly precipitation [LT−1]; β e is the evaporation rate of the groundwater [dimensionless]; Q w n is the monthly actual evapotranspiration [LT−1]; Q c n is additional monthly groundwater recharge and discharge from such activities as local pumping [LT−1]; S y e is the comprehensive specific yield [dimensionless]; and ∆ t is a time step chosen to be 1 month, which is the same as the time interval of the GRACE data [T].
The GWSA ( G W g ) acquired from GRACE indicates the difference between the groundwater storage in a certain month and the average groundwater storage from January 2004 to December 2009, and the value is presumed to be correct. The equation between G W g and groundwater level (GWL) is as follows:
G W g e n = S y e ( h e n − h a v g e )
where h a v g e is the average GWL for cell e from January 2004 to December 2009 (L).
Combining Equation (2) and Equation (3), Equation (1) can be adapted to the following:
G W g e n + 1 ( − ∑ i T e i L e i S y e ( d e i 1 + d e i 2 ) − A e ∆ t ) = − ∑ i G W g e n + 1 T e i L e i S y e i ( d e i 1 + d e i 2 ) − ∑ i T e i L e i ( h a v g e i − h a v g e ) ( d e i 1 + d e i 2 ) − A e ∑ e Q n − G W g e n ∆ t A e
Due to the restriction of in situ wells that are nonuniformly distributed in a majority of areas and the difficulty in achieving the initial hydraulic gradient in groundwater systems, it is assumed that the initial hydraulic gradient can be estimated by multiplying the surface slope by an unknown coefficient, which is expressed as the hydraulic gradient coefficient:
T e i h a v g e i − h a v g e d e i 1 + d e i 2 L e i ≈ T e i Z a v g e i − Z a v g e d e i 1 + d e i 2 C h L e i
where Z a v g e and Z a v g e i are the average surface elevations at cell e and the neighboring cell e i , respectively, and C h is the hydraulic gradient coefficient.
The boundary conditions of the groundwater flow model are set as Dirichlet-type boundaries, which means that the groundwater storage changes at the boundaries are known. Then, the G W g for each cell can be acquired from the equations that can be resolved by considering the boundaries and initial conditions.

3.2. Methods for Bridging the Gap between GRACE and GRACE-FO Data

Numerous studies have extensively examined the assessment of groundwater storage variations in the NCP utilizing GRACE data. Nonetheless, there has been minimal discourse regarding the data gap during the transition from the GRACE to GRACE-FO satellites, spanning from July 2017 to May 2018. The absence of data during this gap period has significantly impeded the ability to monitor the uninterrupted trend of groundwater storage changes. Li et al. [32] compared different driven data methods and found that the approach composed of principal component analysis (PCA), least squares (LS), and multiple linear regression was the most robust alternative to predicting and reconstructing the grid-based TWSC for global watersheds, and its reconstructed TWSC correctly reproduced the EI Niño-Southern Oscillation signal; however, its disadvantage was that it applied only to 26 watersheds. Therefore, Li et al. [33] established a global-scale long-term (1979–2020) TWSA (spatial scale: 0.05°) in 2021. Moreover, Mo et al. [34] devised a novel approach utilizing a Bayesian convolutional neural network to reconstruct missing signals based on hydrological and climatic inputs. Subsequently, they employed the reconstructed data alongside existing signals to investigate regions experiencing persistent water scarcity during the period of 2017 to 2018. Concurrently, data with a 1° × 1° spatial resolution were used in this study for TWSA data comparison.
In this study, we employed the calibrated groundwater storage model to simulate groundwater storage changes during the window period. However, we needed to apply other data sources of groundwater storage changes to carry out the comparison and determine the boundary conditions during the window period. Data standardization encompasses the transformation of data into a universally adopted format or structure, guaranteeing coherence and uniformity across diverse data origins. It entails the establishment of a comprehensive set of regulations, principles, and protocols governing the collection, organization, and representation of data. The significance of standardization in data management and analysis lies in its ability to enhance data quality, promote interoperability, and enable proficient decision-making processes. When data maintain consistent structuring and formatting, they facilitate the identification of patterns, trends, and interrelationships among datasets. In this study, we employed Equation (6) to standardize Li et al.’s data, aiming to harmonize the reconstructed data with the GRACE-retrieved data. This standardization process resulted in a reduction in the disparity between the GRACE-retrieved datasets, enabling their incorporation into the model.
s c a l e d _ d a t a i j = t a r g e t _ d a t a i j − t a r g e t _ m e a n t a r g e t _ s t d × r e f _ s t d + r e f _ m e a n   ( i = 2003 ,   … ,   2020 ; j = 1 , … , 12 )
where s c a l e d _ d a t a i j refers to the data after standardization from year i and month j , t a r g e t _ d a t a i j is the data before standardization from year i and month j , t a r g e t _ m e a n refers to the mean value of the target data, t a r g e t _ s t d is the standard deviation of the target data, r e f _ s t d is the standard deviation of the reference data, and r e f _ m e a n is the mean value of the reference data.

3.3. GRACE Groundwater Storage Drought Index (GGDI)

Primarily, a monthly climatology for C i was determined:
C i = 1 n i ∑ 1 n i ( G W S A i ) ( i = 1 ,   2 ,   … ,   12 )
where C i is the climatology for month i , and n refers to the number of years from 2003 to 2020.
In this study, the application of monthly climatology enabled the separation of the influence of seasonality on groundwater storage changes [35]. Next, by subtracting monthly climatology from GWSA, the groundwater storage deviation (GSD) was derived as the overall deviation in the volume of groundwater storage, accounting for seasonal variations. Ultimately, the groundwater storage deviation (GSD) was normalized by subtracting the mean and dividing it by the standard deviation, thereby standardizing the GSD values:
G G D I = G S D t − x ¯ G S D s G S D
where x ¯ G S D is the mean value, and s G S D refers to the standard deviation.
The GGDI signifies the normalized net deviation in the volume of groundwater storage, and it was employed to assess the characteristics of groundwater drought. The current study utilized the nondimensional GGDI to examine groundwater drought attributes. The evaluation criteria for GGDI are shown in Table 1 [36].

3.4. Model Evaluation

In this study, the performance of the model was evaluated using an objective function known as the total root mean square error ( ϕ ). The total RMSE serves as a measure of the overall agreement between the observed and computed GWSA values. This objective function compares the observation values ( G W g _ o b s ) of GWSA retrieved by GRACE with the computed values ( G W g _ s i m ) of GWSA obtained from the model:
ϕ = 1 n ∑ m = 1 N ∑ n = 1 N d t ( G W m g _ s i m n − G W m g _ o b s n ) 2
where N is the number of cells with an unknown GWSA in the model area. The time step is N d t , and m and n are the continuous counts of grid cells and time steps, respectively.
An optimization algorithm of the compound optimization strategy called the shuffled complex evolution algorithm (SCE-UA), which was proposed by Duan et al. [37], was used to minimize ϕ . The SCE-UA algorithm combines the advantages of the genetic algorithm and simplex method, and it can quickly converge to the global optimal solution by combining deterministic compound search technology with the principle of biological competitive evolution in nature.
In this study, three methods were employed to assess the deviation between simulated groundwater storage values and observed values: correlation coefficient (CC) analysis and root mean square error (RMSE) and Nash-Sutcliffe efficiency coefficient (NSE).
C C = ∑ n = 1 N d t ( X n − X ¯ ) ( Y n − Y ¯ ) ∑ n = 1 N d t ( X n − X ¯ ) 2 ∑ n = 1 N d t ( Y n − Y ¯ ) 2
R M S E = ∑ n = 1 N d t ( X n − Y n ) 2 N d t
N S E = ∑ t = 1 T ( X n t − Y n t ) 2 ∑ t = 1 T ( X n t − Y ¯ ) 2
where X n is the simulated value, Y n is the observed value, X ¯ is the mean of the simulated values, X n t refers to the observed value at time t, Y n t is the simulated value at time t, and Y ¯ is the mean of the observed values.

3.5. Dataset Preparation

The data used in this study included terrestrial water storage (TWS), soil moisture (SM), snow water equivalent (SWE), precipitation data from three diverse datasets, actual evapotranspiration (AET) data from three different datasets, and in situ observation data (Table 2). There was a period of missing data due to the transition from the GRACE to GRACE-FO satellites; therefore, the TWS data from July 2017 to May 2018 were missing. However, other short-range missing data were interpolated by the simple average of the two months before and after the month [38].

3.5.1. Precipitation and Evapotranspiration Data

In this study, three precipitation datasets were utilized: CMA data obtained from the China Meteorological Administration, Tropical Rainfall Measuring Mission 3B43 (TRMM) data, ERA data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF), and Peng datasets (PENG) [39]. The time step for these datasets was monthly. Additionally, three datasets of actual evapotranspiration were employed: MOD16 data from the Moderate Resolution Imaging Spectroradiometer (MODIS), ERA5 data, and Global Land Evaporation Amsterdam Model v.3.5a (GLEAM v.3.5a) data. The time step for these AET datasets was also monthly.
Figure 2 presents a comparison of different precipitation and AET datasets. Among the various precipitation datasets, ERA5 had higher values than the others. In terms of AET data, MOD16 AET data were higher, while GLEAM AET data were smaller. By examining and comparing these datasets, for the subsequent discussions, the study primarily relied on TRMM precipitation data and ERA5 AET data for further analysis and discussion.

3.5.2. GRACE-Derived Data

The present study considered the GRACE monthly 0.5° Level-3 mascon data (RL 06) processed by the Jet Propulsion Laboratory (JPL). The time period applied was from January 2003 to December 2020. The SWE and SM were obtained from Global Land Data Assimilation System (GLDAS) Noah V2.1. To ensure consistency with the GRACE data and establish a comparable interpretation, the soil moisture and snow water equivalent data were transformed into soil moisture anomaly (SMA) and snow water equivalent anomaly (SWEA) by subtracting the average values from the period of 2004 to 2009. This approach allowed for a meaningful comparison between the anomalies derived from the GRACE data and the anomalies calculated from the soil moisture and snow water equivalent datasets. Given that in the NCP, there is a dearth of surface water in arid regions, we opted to disregard the surface water storage capacity of rivers, lakes, and reservoirs [24]. Consequently, the GWSA was determined by subtracting the soil water anomalies and snow water equivalent anomalies within the region from the TWS anomalies. All the data used to construct the model were upscaled to a spatial resolution of 1° in the beginning of the study.
The backward difference in the GWSA in the following formula was used to calculate the GWS changes [23]:
Δ G W S Δ t = G W S A ( t ) − G W S A ( t − 1 ) Δ t
where t refers to the sequence in the monthly GWSA series, and Δ t is evaluated as one month to be in accordance with the temporal resolution of the GRACE observations.

3.5.3. In Situ Data

For the purpose of validating the results obtained from the model calculations, daily in situ groundwater level data (in meters) located in the Beijing-Tianjin-Hebei region were collected [24]. Wells with relatively complete and reliable data were carefully selected for inclusion in the study. To ensure consistency, the data were resampled on a monthly basis, allowing the use of authentic well data for the analysis and evaluation of the model’s performance. To ensure consistency in data interpretation, the groundwater level data were transformed into groundwater level anomalies (GWLA) by subtracting the average value from the period of 2004–2009. This approach aligned the GWLA data with the same contextual meaning as the GWSA data.

3.6. Model Development

Since the GWSA retrieved by GRACE and GLDAS represents groundwater storage anomalies of the entire aquifer thickness, the model is two-dimensional, which means, vertically, the model considers only a one-layer structure. The range of the model was 112°E to 121°E and 33°N to 42°N, and the area in this study was discretized into 81 grid cells with a size of 1°. The external grids of the model were considered as the Dirichlet boundary (cells G1−G10, G18−G19, G26−G28, G34−G37, G44−G46, G54−G55, G63−G64, G72−G81), and the GWSA in 45 internal grid cells was simulated (Figure 3a). In this study, the main recharge source was precipitation infiltration, and the main discharge items were evapotranspiration and groundwater exploitation. The changes in surface water were quite small and were neglected. More than 70% of groundwater was used for irrigation; therefore, a substantial part of irrigation water was lost by ET. Groundwater discharge can be determined by the actual evapotranspiration and the coefficient [23].
The model parameters contained the transmissivity ( T ), comprehensive specific yield ( S y ), precipitation infiltration coefficient ( α ), evapotranspiration coefficient ( β ), and hydraulic gradient coefficient ( C h ). The research area was generalized into three zones on the basis of the hydrogeological characteristics: a high-permeability alluvial-proluvial pore aquifer in the accumulation plain; a medium-permeability alluvial-proluvial pore aquifer in the accumulation plain and coastal plain; and a weak-permeability alluvial-proluvial pore aquifer in the accumulation plain and coastal plain (Figure 3a). Figure 3b describes a hydrogeological cross section from A to A’ (the geographic location is shown in Figure 3a) in the NCP [25]. In addition, the hydraulic parameters of the Quaternary aquifer in the NCP and shallow groundwater flow field in natural state in the NCP are shown in Figures S1 and S2 [40]. Equation (14) can be used to estimate the parameters for each grid when obtaining the parameter values for each zone.
P G i = ∑ k = 1 3 A r e a k × P k
where P G i refers to the grid G i parameter, A r e a k refers to the area percentage for different aquifers, and P k refers to the parameters (transmissivity, comprehensive specific yield, precipitation infiltration, and evapotranspiration coefficient) for the k th type of aquifer.

4. Results

4.1. Hydrogeology Parameter Estimation Using the Groundwater Storage Model

In the present study, the model utilized the GWSA from the GRACE and GLDAS inversion in January 2003 as the initial condition. The time-series GWSA of the GRACE and GLDAS inversion in the Dirichlet boundary grids (Figure 3a) was used as the data boundary. The calibration period spanned from January 2003 to June 2017, while the validation period spanned from June 2018 to December 2019, with a monthly time step. The study employed Python programs based on the SCE-UA algorithm to estimate the optimal parameters by setting parameter thresholds. Table 3 provides the CC, RMSE and NSE, while Figure 4 illustrates the comparison between the GWSA obtained from the GRACE and GLDAS inversion and the simulated GWSA. During the calibration period, the correlation coefficients of the grid cells within the NCP region generally exceeded 0.85. In the validation period, the correlation coefficients remained mostly above 0.4. However, the RMSE in all grid cells was higher than the 2 cm equivalent water height (EWH), possibly due to the influences of coastal area data and boundary conditions. In general, the observed trend during the validation period aligned well with the trend during the simulation period.
The present study executed a sensitivity analysis using the Morris method [24] to assess the impact of parameter uncertainty on the model. The Morris method involves systematically varying the input parameter values of the model and examining the resulting changes in model outputs. In this case, 1000 groups of parameter samples were generated, and the mean values and standard deviations were used to create a logarithmic coordinate chart to visualize the results of the Morris analysis (Figure S3). The sensitivity analysis revealed that the hydraulic gradient coefficient had the highest sensitivity, indicating that it had a significant influence on the model outputs. The hydraulic gradient coefficient represents the state of groundwater flow patterns. If the distribution of the average groundwater level is estimated, the hydraulic gradient coefficient can be calculated by using Equation (5). However, in most cases, the true distribution of groundwater level remains unknown. The hydraulic gradient coefficient as exhibiting the highest sensitivity shows groundwater flow patterns importantly affect the model accuracy. The specific yield showed slightly lower sensitivity than the hydraulic gradient coefficient, while the hydraulic conductivity had the lowest sensitivity among the parameters examined.
Furthermore, the study investigated the influence of multisource remote sensing data products on the model simulation results. Various combinations of precipitation and evapotranspiration remote sensing data were used as model inputs. It is worth noting that the optimized hydrogeological parameters remained unchanged during this analysis. The simulation results from each data combination were compared, and it was observed that the variation trends of the results were remarkably consistent across all groups. Although the RMSE values were slightly higher for some combinations, the correlation coefficients were very similar and did not exhibit any notable differences (Figure S4). This result suggests that while there may be numerical discrepancies among different remote sensing data, the calculated GWSA was not significantly affected by the choice of remote sensing data combination.
Sun et al. [24] introduced the notion of comprehensive specific yield, and their findings strongly supported the notion of an optimized comprehensive specific yield, which appears to be reasonable. This article employed the concept of the comprehensive specific yield. Figure S5 shows the spatial distribution of relevant hydrogeological parameters in the NCP after downscaling. The minimum value of the precipitation infiltration coefficient was 0.016, and the maximum value was 0.323. Its distribution range gradually decreased from the piedmont plain to the coastal plain. The precipitation infiltration coefficient in the piedmont plain was the highest, ranging from 0.255 to 0.323. The range of the central and eastern plains was 0.135–0.255, with the lowest range of 0.050–0.135 in the coastal plain and 0.015–0.050 in a few areas. The minimum value of the evaporation coefficient was 0.07, and the maximum value was 0.661. The evaporation coefficient of the piedmont plain was the highest, ranging from 0.500 to 0.661. The range of the central and eastern plains was 0.050–0.250, with a few areas ranging from 0.250 to 0.500. The range of comprehensive specific yield was 7 × 10−5−0.0012, and the distribution range gradually decreased from the piedmont plain to the coastal plain. The obtained results exhibited a high degree of similarity and comparability to those reported in previous studies [13,24,41].

4.2. Downscaled GWSA Changes

In the downscaling preparation of the model, the outermost grid cells of the study area remained unchanged, preserving the original boundaries and initial conditions derived from the GWSA retrieved by GRACE and GLDAS. The inner region of the study area, consisting of 45 grid cells, underwent further subdivision into a higher resolution grid system. Specifically, these 45 grid cells were subdivided into a total of 19,600 smaller grid cells. This refinement ranged from 1° grid cells (Figure 5a) to more detailed 0.05° grid cells (Figure 5c). The precipitation and actual evaporation datasets in 0.05° model were the PENG dataset and MOD16 dataset, respectively. Both datasets had a spatial resolution of 0.05°, matching the resolution of the downscaling model. The hydrogeological parameter values used in the downscaling model were processed using Equation (14) with the same spatial resolution as the downscaling model. The downscaling model was driven based on the precipitation and actual evaporation data with a spatial resolution of 0.05°.
The downscaled model simulation period spanned January 2003 to December 2020, regardless of the time gap between the GRACE and GRACE-FO satellites. It was evident that the general spatial distribution trend of the downscaled GWSA (Figure 5d) closely resembled the GWSA retrieved by GRACE and GLDAS (Figure 5b). In addition, the downscaled GWSA exhibited finer spatial characteristics, with smoother and more detailed transitions between regions. By examining the time-series trend of the downscaled results in Figure 6, it was observed that the change trend of the downscaled GWSA in the grid cells within the NCP region was consistent with that of the GWSA derived from GRACE and GLDAS. This result indicates that the downscaled model effectively captured and reproduced the temporal variations in groundwater storage in the NCP region, aligning with the findings obtained through the GRACE and GLDAS inversion.
In this study, a total of 12 wells with complete and reliable data were selected for comparison with the downscaled GWSA. Figure 7 presents the comparison results between the downscaled GWSA and the observed GWLA data. Table 4 shows the correlation coefficients between the downscaled GWSA and the observed GWLA data. The results demonstrate a close alignment between the trends of the downscaled GWSA and the GWLA data. Furthermore, the correlation coefficients between the variation trends of the downscaled GWSA and GWLA were consistently greater than 0.5, with the highest correlation coefficient reaching 0.72, which indicates that the downscaled GWSA effectively reflects the dynamic changes in groundwater within a small range in the NCP. Overall, the downscaled GWSA serves as a reliable indicator of groundwater dynamics in the study area.

4.3. Bridging the Gap between GRACE and GRACE-FO Data

The TWSA data for the NCP from January 2003 to December 2020 were extracted from the studies conducted by Li et al. [33] (Li), Mo et al. [34] (Mo), and the CLSM model (CLSM), and the data were compared with the TWSA retrieved by GRACE. Note that the TWSA of Li spanned from January 2003 to June 2020, and the TWSA of Mo and the GLDAS CLSM model spanned from January 2003 to December 2020. As shown in Figure 8, Li data have the same trend as GRACE TWSA data, and its correlation coefficient reached 0.96, while Mo and CLSM TWSA data have shown opposite trends to GRACE TWSA data since 2008. In fact, after 2008, contrasting trends emerged, and the mean numerical disparity between the Li dataset and GRACE dataset at a given time point amounted to 8 cm EWH.
The average difference between the Li standardized TWSA data and the GRACE-retrieved TWSA data was reduced to 2 cm EWH. The snow water equivalent and soil moisture were subtracted from the standardized TWSA to obtain the GWSA. The GWSA from July 2017 to May 2018 was selected as the input for the model, while the GWSA for the remaining time periods still used GRACE-retrieved data. The calculated results are shown in Figure 9. Compared with the calculated GWSA of Mo and CLSM, Li data and GRACE data have better trends and higher correlations. Given the groundwater storage changes in the model boundaries, the groundwater storage changes in the NCP for the gap period can be calculated by using numerical groundwater storage model. In this study, Li data (July 2017–May 2018) is referred to produce the model boundary information.

5. Discussion

5.1. GWS Changes in the NCP and Five Administrative Regions

The NCP covers parts of Beijing (BJ), Tianjin (TJ), Hebei (HB), Henan (HN), and Shandong (SD). Figure 10 depicts the time series diagram of the monthly GWSA data and precipitation data of the NCP and its five subregions from 2003 to 2020, and Table 5 shows a summary of the comparison of the results for different time periods and research regions. From 2003 to 2010, the rate of decline of the GWSA in the NCP was 0.91 cm/yr, which was lower than that presented in the research results of Feng et al. [5]. The possible reason is that the study area of Feng et al. was 370,000 km2, which was much larger than the study area of this study. The loss rate of GWSA from 2003 to 2014 in the NCP was 1.40 cm/yr. The result from Feng et al. [15] from 2002 to 2014 obtained through GRACE was 7.2 ± 1.1 km3/yr, which was still greater than the result in the present study after being converted to the same units. One potential explanation is the substantial extent of the region, spanning a large area of 320,000 km2, which surpassed the scope of the current study. While the research findings displayed slight numerical discrepancies, the time series of GWSA remained fundamentally congruent, depicting a discernible seasonal variation pattern. In the same research period, the result from Liu et al. [42], whose study was conducted within an identical study area as the present study, reported a groundwater depletion rate of 1.66 ± 0.17 cm/yr, demonstrating a high level of similarity. Furthermore, both studies identified primarily similar trends in the variation in groundwater levels. In addition, in Liu et al. [42], the R2 value of the monthly GWSA from 2005 to 2016 and the in situ GWSA derived from GRACE was 0.91, and the variation trend (−2.21 ± 0.15 cm/yr) exhibited a remarkable resemblance to the result (−1.91 cm/yr) obtained in the current study. In this research, the rate of the decline of the GWSA from 2003 to 2015 in the NCP was 15.2 mm/yr. From 2005 to 2013, the rates of decline of the GWSA of the HB and TJ were 15.9 mm/yr and 18.6 mm/yr, respectively. These results were fundamentally similar to the results obtained by Gong et al. [14]. Xu et al. [43] also obtained GWSA data from GRACE from 2003 to 2017, and the correlation coefficient between the GRACE-retrieved GWSA (−19.96 ± 3.6 mm/yr) and the field-measured data (−16.80 ± 1.00 mm/yr) reached 0.75. In this research, the alteration rate of the GWSA during the identical time period was −17.15 mm/yr, which was remarkably consistent. In the research on groundwater depletion in the NCP by Zhao et al. [44], the depletion rate of groundwater storage from 2004 to mid-2016 was 1.7 ± 0.1 cm/yr, and from mid-2013 to mid-2016, the depletion rate was 3.8 ± 0.1 cm/yr. The rates of decline in the present study from 2004 to mid-2016 and from mid-2013 to mid-2016 were 1.76 cm/yr and 2.47 cm/yr, respectively. The value from the middle of 2013 to the middle of 2016 was slightly smaller, but they all revealed that the groundwater decline in the NCP was more severe after 2013, reflecting the long-term decline in groundwater storage. From 2003 to 2016, the variation rate of the GWSA in the NCP was −16.1 mm/yr, which was quite close to the result (−17.2 ± 0.8 mm/yr) of Zheng et al. [30]. In general, the time series of downscaled GWSA effectively captures the fluctuations in groundwater storage within the NCP as well as in certain administrative regions.
As shown in Figure 10, the GWSA in the NCP and five subregions presented a continuous downward trend. Moreover, the time series of the GWSA showed obvious seasonal fluctuations. The GWSA value declined rapidly in spring and summer and increased slowly in autumn and winter. The NCP is an area with two types of crops per year. One is winter wheat, which is sown in October and harvested in June of the next year, and the other is summer maize, which is sown in July and harvested in October of the same year [45]. The seasonal fluctuations in the GWSA were mainly caused by a large amount of agricultural irrigation pumping groundwater and precipitation [42]. The growth of wheat depends heavily on groundwater irrigation, while corn requires little groundwater because of the intensive precipitation during the summertime [31,46]. There was also a time lag between changes in the downscaled GWSA and precipitation. This indicated the process of groundwater recharge by precipitation [47]. The GWSA decreased most when the precipitation was highest, which implies that precipitation is not the main influencing factor of groundwater storage changes. Zhuo et al. [48] declared that non-climatic factors play a leading role in the change in shallow groundwater storage in the NCP. Leng et al. [49] also reported that groundwater irrigation is the main reason for the decline in groundwater level in the main agricultural areas in the NCP, and climate change has little impact on the change in groundwater level.
Table 6 shows the four-stage change rate of the downscaled GWSA in the NCP and five administrative regions. During the period spanning from January 2003 to December 2008, the GWSA within the five administrative regions exhibited a modest declining trend. However, from December 2008 to June 2009 and from October 2010 to June 2011, the rate of decline of the GWSA increased significantly (Figure 10). This may be caused by two extreme drought events in the North China Plain [13]. Generally, meteorological drought is likely to be the main reason for the decrease in soil water and the surface water resource deficiency [50], which will increase the demand for agricultural water, leading to more intensive groundwater extraction [14,31]. Between January 2009 and December 2014, there was a notable acceleration in the rate of GWSA decline. Particularly noteworthy is the fact that within this period, the HN region experienced a remarkable surge in the rate of decline, reaching 32 times higher than the previous cycle. In 2013 and 2014, the precipitation in the NCP was less than the average of previous years (549 mm), with precipitation of 521 mm and 469 mm, respectively. Drought occurred from August 2013 to September 2014, and the most severe drought occurred in July 2014 [36]. In the wake of the intense drought, the rate of decline of the GWSA experienced a swift escalation within the NCP. During the period from August 2014 to May 2015, the GWSA in the region decreased by 54 mm. From January 2015 to December 2017, the rate of decline of the GWSA reached the maximum of the four stages. The loss rate of GWSA in the SD district reached 0.245 cm/month. The major cause may be that BJ and TJ are the primary cities that receive replenishment from the SNWTP [50]. Critical drought occurred in HN, SD and other areas after 2014, resulting in the massive exploitation of groundwater, which led to a serious decline in groundwater storage [50,51]. From January 2018 to December 2020, the overall rate of decline of the GWSA within the NCP exhibited a marked deceleration, with a decrease of nearly half. Especially in BJ and TJ, the variation trend of the GWSA became positive, with values of 0.095 cm/month and 0.048 cm/month, respectively. The SNWTP and other management measures used to control groundwater overexploitation have been effective. The SNWTP was officially opened to water in December 2014. It aims to alleviate the shortage of water resources in the central and northern regions of China. The water transfer plan assumes that from 2010 to 2025, the urban domestic and industrial water demand will remain unchanged, and the main goal is to reduce the water consumption of groundwater in urban domestic and industrial areas [29]. Nevertheless, with the burgeoning development of urbanization and industry, the demand for urban domestic and industrial water is nondeterministic [52]. This means the amount of surface water redistributed to agricultural irrigation is still as uncertain as before [29]. Water diversion and precipitation alone cannot effectively solve the decline in groundwater storage, and it will take a considerable amount of time to recover from the groundwater overexploitation that occurred over the past few decades [47]. In addition, the auxiliary projects of the SNWTP were still in the stage of gradual completion from 2013 to 2017, resulting in its actual water supply volume being lower than the average designed water supply volume for many years at the initial stage of the project [41].

5.2. Spatial Variation in GWS Changes

Compared to the results of previous research, the time series of the downscaled GWSA effectively captured the changes in groundwater storage within the NCP and its five administrative regions. Figure S6 demonstrates the spatial distribution of the monthly average GWSA in the NCP in two periods (from 2003 to 2020). The GWSA presented a seasonal trend. In the fall and winter seasons, the GWSA variation range was comparatively small, while in the spring and summer seasons, the GWSA variation was relatively large, which was very similar to the results of Xu et al. [43] and Liu et al. [42]. Moreover, the spatial decline trend of groundwater storage gradually increased from northeast to southwest, which was similar to the result of Xu et al. [43].
Figure 11 shows the spatial distribution of the change rate of the GWSA in four stages in the NCP. From December 2008 to January 2009, severe drought events occurred in Northern China [13]. In December 2014, the Mid-Route of the SNWTP officially supplied water to Beijing, Tianjin, Hebei, and Henan. Therefore, the present study took 2009 and 2015 as time nodes and calculated the monthly groundwater storage change rates from 2003 to 2008, 2009 to 2014, January 2015 to June 2017, and June 2018 to December 2020 using the linear regression method. The spatial distribution of groundwater storage change trends in the NCP is shown in Figure 11a–d. From 2003 to 2008, the groundwater storage in the NCP, except the southern region, showed a slight decline. The regions with heavy downward rates were distributed mainly in Tianjin, Cangzhou, and Xingtai. From 2009 to 2014, the groundwater storage in the NCP showed an obvious downward trend compared with the previous period, and the rate of decline gradually increased from the northeast to southwest. The total amount of groundwater exploitation remained high without decreasing; therefore, the water level continued to decline at an accelerated rate. Xu et al. [43] found that from 2003 to 2017, the groundwater storage in the NCP showed a long-term downward trend, and the downward trend gradually increased from the northeast to southwest, which was similar to the results in this study. From January 2015 to December 2017, the rate of decline of groundwater storage in the NCP continued to increase. The areas with the most serious decline were distributed mainly in Xingtai city and Handan city. In the study conducted by Zhao et al. [44], the rate of groundwater storage depletion in the NCP intensified after 2013. Furthermore, through the analysis of vertical displacement data from 116 sediment GPS stations, it was discovered that land subsidence experienced an accelerated pace after 2013, which aligned closely with the increased rate of groundwater depletion. From January 2018 to December 2020, the areas with relatively large groundwater decline rates in the NCP shifted from the northeast and southeast regions to the southwest regions, especially in Xingtai, Handan, and Henan. The spatial distribution was highly similar to the results of Sun et al. [24]. The groundwater depletion in the NCP achieved obvious results after treatment. In particular, groundwater in the Beijing area is recovering [53]. Numerous measures, including the regulation of groundwater overexploitation, water supply from the SNWTP, ecological augmentation of rivers and lakes, and agricultural water conservation, have yielded initial outcomes. However, the magnitude of groundwater drawdown cones, formed due to the prolonged imbalance between groundwater exploitation and recharge, remained substantial. The average groundwater level of shallow aquifers in the NCP declined by more than 15 m and that of deep aquifers declined by more than 28 m over the past 50 years because of long-term and continuous groundwater overexploitation [29]. Moreover, as the supporting project of the SNWTP was still in the stage of gradual completion, its actual water supply was lower than the average expected water supply for many years; therefore, it was difficult to alleviate the serious problem of groundwater depletion in a short time. Nevertheless, as the SNWTP continues, the actual water supply from the project is expected to increase steadily over time. Its impact on promoting groundwater recharge in the NCP is anticipated to become increasingly evident in the future [41].

5.3. Estimation of the GGDI from a Spatiotemporal Perspective

Thomas et al. [35] proposed a method to assess the incidence of groundwater droughts using data from the GRACE satellite mission. They demonstrated that the normalized groundwater storage deviations derived from GRACE could be used to measure groundwater storage deficits during the GRACE era, and they called this the GGDI. Numerous researchers have conducted extensive investigations and analyses on the spatial distribution, temporal progression, and trend characteristics of drought in the study area using the GGDI. The findings consistently demonstrate the reliability and robustness of GRACE-derived quantitative outcomes for drought assessment [36]. Consequently, this study utilized the downscaled GWSA to compute the GGDI, enabling an exploration of the spatiotemporal trends associated with the GGDI.
Figure 12 depicts the spatial distribution of the GGDI in the NCP. Table 7 presents the rate of GGDI change across four distinct periods. From 2003 to 2008, with the exception of Henan, the GGDI in the remaining three regions exhibited a marginal decline, and the trend was basically consistent. Between 2009 and 2014, there was a notable increase in the rate of decline, as observed in the GGDI time series. Notably, the Henan region experienced the highest GGDI variation rate, at −0.0192. From a spatial perspective, the most pronounced decrease in the GGDI was observed in the southern areas of the North China Plain, specifically in Henan, Handan, and Liaocheng. During the period of 2015 to 2017, the NCP experienced the highest rate of decline in the GGDI, reaching 0.0206. Notably, Beijing exhibited the highest rate of decline among the five administrative regions, reaching 0.0301, which was almost 0.1 higher than the monthly average rate of decline of the GGDI in the NCP. In contrast, the rate of decline of the GGDI in the Henan region was relatively lower during this period. From a spatial perspective, the pronounced rate of decline of the GGDI was primarily concentrated in the central region of the North China Plain, with Beijing experiencing a particularly severe decline. Between 2018 and 2020, the rate of decline of the GGDI in the NCP dropped by almost half. Notably, the GGDI change rate in the Beijing and Tianjin regions may even exhibit positive values, indicating a potential improvement in groundwater drought. However, the rate of decline of the GGDI in the Henan region continues to increase, reaching a maximum value of −0.0274 during the four stages. From a spatial perspective, the areas with the most prominent decline in the GGDI within the NCP were situated in the Piedmont plain of Hebei and Henan. The spatial distribution of the GGDI change rate in the four stages was very similar to the spatial distribution of the GWSA change rate in Figure 12. Since 2015, the NCP has experienced a period of mild groundwater drought, which escalated to a moderate level after 2016. Despite the opening of the middle line of the SNWTP for water supply in December 2014, it has not rapidly alleviated the ongoing groundwater drought trend in the NCP, mainly due to the water supply not meeting the anticipated levels [41]. The spatial distribution pattern of the GGDI changes from 2018 to 2020 suggested an improvement in the groundwater drought situation in the Beijing-Tianjin area, possibly due to the influence of the SNWTP and other water-saving measures. However, it is important to note that the GGDI may exhibit a delayed response to implemented mitigation measures. Previous studies have also mentioned similar conclusions [35].

5.4. Limitations and Perspectives

The utilization of GRACE data has been widely validated in numerous studies as an effective means to investigate regional groundwater storage changes. It offers a swift and cost-free approach for monitoring such changes. However, it is important to acknowledge that estimates of GWS changes derived from GRACE data still entail a degree of uncertainty. This uncertainty arises from errors present in the TWS changes derived from GRACE and in other components obtained from the GLDAS hydrological model. In this study, the components utilized were derived from the GLDAS hydrological model. While this model may not be able to capture the short-term impacts of human activities on groundwater storage changes, it reliably reflects long-term trends [23].
In this study, the analysis focused solely on precipitation and actual evapotranspiration as the source and sink components, respectively, without adequately considering additional factors such as the SNWTP and groundwater extraction. It is crucial to note that the NCP, being a densely populated region in China, experiences significant influences on groundwater storage changes due to various human activities. Furthermore, the model employed in the study used relatively coarse parameter divisions, relying primarily on the topography and permeability coefficient distribution within the study area. Given the complexity of groundwater dynamics and the numerous anthropogenic factors at play, a more comprehensive assessment incorporating a wider range of factors and more refined parameter divisions would provide a more accurate representation of the groundwater storage changes in the NCP. However, the land subsidence in the NCP caused by years of groundwater overexploitation has caused changes in the hydrogeological conditions in some areas [14]. After downscaling the GRACE gap period data of Li et al. [33] as a boundary condition, good results were obtained; however, the lack of in situ observation data made it difficult to verify the reliability of the downscaled data. If more of the data required for numerical simulation can be obtained in subsequent research, the model can be further improved, and the research content can be further studied.

6. Conclusions

In this study, a groundwater storage change model was constructed in the NCP using a dynamic downscaling model with the GWS inversed from GRACE and GLDAS data as the fitting target, and the SCE-UA algorithm was used to optimize the hydrogeological parameters. The optimized hydrogeological parameters were used to improve the resolution of GWSA data from 1° to 0.05°, and the reliability of the data was verified. The data gap that occurred between the GRACE and GRACE-FO satellites was bridged. By comparison with previous research results, it was found that the model data had high reliability after downscaling. The conclusions are as follows:
(1)
The established groundwater storage model using multiple remote-sensing data demonstrated perfect performance after model calibration and verification. The correlation coefficients between the simulated and GRACE-observed GWSA in the calibration period were all greater than 0.85, and 55% of the correlation coefficients in the validation period were greater than 0.50. The uncertainty analysis of the model showed that the combinations of precipitation and actual evapotranspiration data from different sources had no significant impact on the simulated GWSA outputs. The sensitivity of the hydraulic gradient coefficient was the highest, while the sensitivity of the specific yield was slightly lower than that of the hydraulic gradient coefficient.
(2)
The downscaled GWSA in the NCP showed a similar and finer spatial distribution when compared with that retrieved by GRACE and GLDAS as well as consistent changes with the in situ observations. Meanwhile, the missing GWSA values during the period of transition between the GRACE and GRACE-FO satellites were bridged. The comparison of the results with previous studies demonstrated favorable performance and was deemed reasonable, affirming the validity and rationality of the model in compensating for the downscaling of data in the empty window period.
(3)
The GWSA changes in the five subregions (BJ, TJ, HB, HN, and SD) showed different patterns from 2003 to 2020. From 2003 to 2008, the GWS fluctuated and declined except in HN. From 2008 to 2014, the GWS declined overall. From January 2014 to June 2017, the GWS showed a rapid downward trend. From June 2018 to December 2020, the downward trend of the GWS was significantly slower than that of the previous stage, and in the BJ region, the variation trend of the GWS showed a slow upward trend. This result may be due to the initial success of the STNWTP and control measures for groundwater overexploitation in the NCP.
(4)
The patterns of the calculated GGDI in the NCP for the time period from 2003 to 2020 were similar to those of the GWSA. The analysis of the GGDI changes in the five administrative regions (BJ, TJ, HB, HN, and SD) over the period from 2003 to 2020 revealed distinct patterns. From 2003 to 2008, the GGDI exhibited fluctuations and an overall decline, except in HN, where it remained relatively stable. Subsequently, from 2008 to 2014, the GGDI showed a general decline across all regions. During the period from January 2014 to June 2017, the GGDI experienced a rapid and significant downward trend. However, from June 2018 to December 2020, the rate of decline of the GGDI slowed notably compared to the previous stage, except in HN. In particular, in the BJ and TJ regions, the GGDI even exhibited a slight upward trend. Overall, the spatial distribution of the GGDI variations closely resembled that of the GWSA.
Overall, the present study provides a reference for the spatial downscaling of GRACE data to achieve a more refined analysis of groundwater storage changes in the NCP.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs15133264/s1, Figure S1: Distribution of hydraulic conductivities of the Quaternary aquifer in the NCP; Figure S2: Shallow groundwater flow field in natural state in the NCP in 1959; Figure S3: Sensitivities of hydrogeology parameters; Figure S4: A comparison of the average GWSA for 12 different scenarios using various driving data in the study area; Figure S5: Spatial distribution of (a) precipitation infiltration coefficients and (b) evapotranspiration coefficient and (c) comprehensive specific yield in the NCP; Figure S6: Monthly average GWSA in the NCP from January 2003 to December 2020 (January-December).

Author Contributions

J.Z.: Writing-original draft, writing-review & editing; L.H.: Conceptualization, methodology, writing-review & editing; J.S.: Data collection & analysis, discussed and writing-review & editing; D.W.: Writing-review & editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant numbers 41877173 and U2167211.

Data Availability Statement

This study provides the model program and input data files as well as the downscaled GWSA datasets, which can be downloaded from https://figshare.com/articles/dataset/Downscaling_GWSA_model_and_datasets/22819976 (uploaded on 15 May 2023) or are available from the first author upon request.

Acknowledgments

We thank the public uncertainty quantification tools from Qingyun Duan at Hohai University and constructive suggestions from three anonymous reviewers. We appreciate the free access of the GRACE mascon solutions and GLDAS, TRMM, ERA5, GLEAM, MODIS, and PENG meteorological and hydrological datasets. GRACE mascon solutions: ftp://podaac-ftp.jpl.nasa.gov/allData/tellus/L3/land%20mass/RL06 (accessed on 11 March 2022); GLDAS: http://disc.gsfc.nasa.gov/hydrology/data holdings (accessed on 11 March 2022); GLEAM: http://www.gleam.eu (accessed on 16 March 2022); MODIS: https://doi.org/10.5067/MODIS/MCD12Q1.006 (accessed on 16 March 2022); TRMM: https://gpm.nasa.gov/category/mission-affiliation/trmm (accessed on 21 March 2022); ERA5: https://cds.climate.copernicus.eu/cdsapp#!/home (accessed on 23 March 2022); PENG: http://poles.tpdc.ac.cn/zh-hans/data/faae7605-a0f2-4d18-b28f-5cee413766a2/ (accessed on 2 April 2022); CMA: http://data.cma.cn/ (accessed on 23 March 2022).

Conflicts of Interest

All authors declare that no conflict of interest exists.

References

  1. The United Nations Educational Scientific and Cultural Organization (UNESCO). The United Nations World Water Development Report 2022: Groundwater: Making the Invisible Visible; UNESCO: Paris, France, 2022; ISBN 978-92-3-100507-7. [Google Scholar]
  2. Masood, A.; Tariq, M.A.U.R.; Hashmi, M.Z.U.R.; Waseem, M.; Sarwar, M.K.; Ali, W.; Farooq, R.; Almazroui, M.; Ng, A.W.M. An Overview of Groundwater Monitoring through Point-to Satellite-Based Techniques. Water 2022, 14, 565. [Google Scholar] [CrossRef]
  3. Dalin, C.; Wada, Y.; Kastner, T.; Puma, M.J. Groundwater Depletion Embedded in International Food Trade. Nature 2017, 543, 700–704. [Google Scholar] [CrossRef] [Green Version]
  4. Chatterjee, R.S.; Pranjal, P.; Jally, S.; Kumar, B.; Dadhwal, V.K.; Srivastav, S.K.; Kumar, D. Potential Groundwater Recharge in North-Western India vs Spaceborne GRACE Gravity Anomaly Based Monsoonal Groundwater Storage Change for Evaluation of Groundwater Potential and Sustainability. Groundw. Sustain. Dev. 2020, 10, 100307. [Google Scholar] [CrossRef]
  5. Feng, W.; Zhong, M.; Lemoine, J.-M.; Biancale, R.; Hsu, H.-T.; Xia, J. Evaluation of Groundwater Depletion in North China Using the Gravity Recovery and Climate Experiment (GRACE) Data and Ground-Based Measurements. Water Resour. Res. 2013, 49, 2110–2118. [Google Scholar] [CrossRef]
  6. López-Alvis, J.; Carrera-Hernández, J.J.; Levresse, G.; Nieto-Samaniego, Á.F. Assessments of Groundwater Depletion Caused by Excessive Extraction through Groundwater Flow Modeling: The Celaya Aquifer in Central Mexico. Environ. Earth Sci. 2019, 78, 482. [Google Scholar] [CrossRef]
  7. Liu, X.; Hu, L.; Sun, K.; Yang, Z.; Sun, J.; Yin, W. Improved Understanding of Groundwater Storage Changes under the Influence of River Basin Governance in Northwestern China Using GRACE Data. Remote Sens. 2021, 13, 2672. [Google Scholar] [CrossRef]
  8. Sheard, B.S.; Heinzel, G.; Danzmann, K.; Shaddock, D.A.; Klipstein, W.M.; Folkner, W.M. Intersatellite Laser Ranging Instrument for the GRACE Follow-on Mission. J. Geod. 2012, 86, 1083–1095. [Google Scholar] [CrossRef]
  9. Rodell, M.; Velicogna, I.; Famiglietti, J.S. Satellite-Based Estimates of Groundwater Depletion in India. Nature 2009, 460, 999–1002. [Google Scholar] [CrossRef] [Green Version]
  10. Hu, L.; Jiao, J.J. Calibration of a Large-Scale Groundwater Flow Model Using GRACE Data: A Case Study in the Qaidam Basin, China. Hydrogeol. J. 2015, 23, 1305–1317. [Google Scholar] [CrossRef]
  11. Tiwari, V.M.; Wahr, J.; Swenson, S. Dwindling Groundwater Resources in Northern India, from Satellite Gravity Observations. Geophys. Res. Lett. 2009, 36, L18401.1–L18401.5. [Google Scholar] [CrossRef] [Green Version]
  12. Li, B.; Rodell, M.; Kumar, S.; Beaudoing, H.K.; Getirana, A.; Zaitchik, B.F.; Goncalves, L.G.; Cossetin, C.; Bhanja, S.; Mukherjee, A.; et al. Global GRACE Data Assimilation for Groundwater and Drought Monitoring: Advances and Challenges. Water Resour. Res. 2019, 55, 7564–7586. [Google Scholar] [CrossRef] [Green Version]
  13. Huang, Z.; Pan, Y.; Gong, H.; Yeh, P.J.-F.; Li, X.; Zhou, D.; Zhao, W. Subregional-Scale Groundwater Depletion Detected by GRACE for Both Shallow and Deep Aquifers in North China Plain. Geophys. Res. Lett. 2015, 42, 1791–1799. [Google Scholar] [CrossRef]
  14. Gong, H.; Pan, Y.; Zheng, L.; Li, X.; Zhu, L.; Zhang, C.; Huang, Z.; Li, Z.; Wang, H.; Zhou, C. Long-Term Groundwater Storage Changes and Land Subsidence Development in the North China Plain (1971–2015). Hydrogeol. J. 2018, 26, 1417–1427. [Google Scholar] [CrossRef] [Green Version]
  15. Feng, W.; Shum, C.; Zhong, M.; Pan, Y. Groundwater Storage Changes in China from Satellite Gravity: An Overview. Remote Sens. 2018, 10, 674. [Google Scholar] [CrossRef] [Green Version]
  16. Zhang, Q.; Li, Y.P.; Huang, G.H.; Wang, H.; Li, Y.F.; Liu, Y.R.; Shen, Z.Y. A Novel Statistical Downscaling Approach for Analyzing Daily Precipitation and Extremes under the Impact of Climate Change: Application to an Arid Region. J. Hydrol. 2022, 615, 128730. [Google Scholar] [CrossRef]
  17. Yin, W.; Hu, L.; Zhang, M.; Wang, J.; Han, S.-C. 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]
  18. Rahaman, M.; Thakur, B.; Kalra, A.; Li, R.; Maheshwari, P. Estimating High-Resolution Groundwater Storage from GRACE: A Random Forest Approach. Environments 2019, 6, 63. [Google Scholar] [CrossRef] [Green Version]
  19. Seyoum, W.; Kwon, D.; Milewski, A. Downscaling GRACE TWSA Data into High-Resolution Groundwater Level Anomaly Using Machine Learning-Based Models in a Glacial Aquifer System. Remote Sens. 2019, 11, 824. [Google Scholar] [CrossRef] [Green Version]
  20. Giorgi, F.; Mearns, L.O. Approaches to the Simulation of Regional Climate Change: A Review. Rev. Geophys. 1991, 29, 191. [Google Scholar] [CrossRef]
  21. Houborg, R.; Rodell, M.; Li, B.; Reichle, R.; Zaitchik, B.F. Drought Indicators Based on Model-Assimilated Gravity Recovery and Climate Experiment (GRACE) Terrestrial Water Storage Observations. Water Resour. Res. 2012, 48, W07525. [Google Scholar] [CrossRef] [Green Version]
  22. Zaitchik, B.F.; Rodell, M.; Reichle, R.H. Assimilation of GRACE Terrestrial Water Storage Data into a Land Surface Model: Results for the Mississippi River Basin. J. Hydrometeorol. 2008, 9, 535–548. [Google Scholar] [CrossRef] [Green Version]
  23. Sun, J.; Hu, L.; Liu, X.; Sun, K. Enhanced Understanding of Groundwater Storage Changes under the Influence of River Basin Governance Using GRACE Data and Downscaling Model. Remote Sens. 2022, 14, 4719. [Google Scholar] [CrossRef]
  24. Sun, J.; Hu, L.; Chen, F.; Sun, K.; Yu, L.; Liu, X. Downscaling Simulation of Groundwater Storage in the Beijing, Tianjin, and Hebei Regions of China Based on GRACE Data. Remote Sens. 2023, 15, 1490. [Google Scholar] [CrossRef]
  25. Cao, G.; Zheng, C.; Scanlon, B.R.; Liu, J.; Li, W. Use of Flow Modeling to Assess Sustainability of Groundwater Resources in the North China Plain. Water Resour. Res. 2013, 49, 159–175. [Google Scholar] [CrossRef]
  26. Shao, J.; Li, L.; Cui, Y.; Zhang, Z. Groundwater Flow Simulation and Its Application in Groundwater Resource Evaluation in the North China Plain, China. Acta Geolo. Sin. Engl. Ed. 2013, 87, 243–253. [Google Scholar] [CrossRef]
  27. Foster, S.; Garduno, H.; Evans, R.; Olson, D.; Tian, Y.; Zhang, W.; Han, Z. Quaternary Aquifer of the North China Plain-assessing and achieving groundwater resource sustainability. Hydrogeol. J. 2004, 12, 81–93. [Google Scholar] [CrossRef]
  28. Ministry of Water Resources of China (MWR). China Water Resources Bulletin 2021. Available online: http://www.mwr.gov.cn/sj/tjgb/szygb/202206/t20220615_1579315.html (accessed on 15 June 2022).
  29. Yao, Y.; Zheng, C.; Andrews, C.; He, X.; Zhang, A.; Liu, J. Integration of Groundwater into China’s South-North Water Transfer Strategy. Sci. Total Environ. 2019, 658, 550–557. [Google Scholar] [CrossRef]
  30. Zheng, L.; Pan, Y.; Gong, H.; Huang, Z.; Zhang, C. Comparing Groundwater Storage Changes in Two Main Grain Producing Areas in China: Implications for Sustainable Agricultural Water Resources Management. Remote Sens. 2020, 12, 2151. [Google Scholar] [CrossRef]
  31. Pan, Y.; Zhang, C.; Gong, H.; Yeh, P.J.-F.; Shen, Y.; Guo, Y.; Huang, Z.; Li, X. Detection of Human-Induced Evapotranspiration Using GRACE Satellite Observations in the Haihe River Basin of China. Geophys. Res. Lett. 2017, 44, 190–199. [Google Scholar] [CrossRef]
  32. Li, F.; Kusche, J.; Rietbroek, R.; Wang, Z.; Forootan, E.; Schulze, K.; Lück, C. Comparison of Data-Driven Techniques to Reconstruct (1992–2002) and Predict (2017–2018) GRACE-like Gridded Total Water Storage Changes Using Climate Inputs. Water Resour. Res. 2020, 56, e2019WR026551. [Google Scholar] [CrossRef] [Green Version]
  33. Li, F.; Kusche, J.; Chao, N.; Wang, Z.; Löcher, A. Long-Term (1979-Present) Total Water Storage Anomalies over the Global Land Derived by Reconstructing GRACE Data. Geophys. Res. Lett. 2021, 48, e2021GL093492. [Google Scholar] [CrossRef]
  34. Mo, S.; Zhong, Y.; Forootan, E.; Shi, X.; Feng, W.; Yin, X.; Wu, J. Hydrological Droughts of 2017–2018 Explained by the Bayesian Reconstruction of GRACE(-FO) Fields. Water Resour. Res. 2022, 58, e2022WR031997. [Google Scholar] [CrossRef]
  35. Thomas, B.F.; Famiglietti, J.S.; Landerer, F.W.; Wiese, D.N.; Molotch, N.P.; Argus, D.F. GRACE Groundwater Drought Index: Evaluation of California Central Valley Groundwater Drought. Remote Sens. Environ. 2017, 198, 384–392. [Google Scholar] [CrossRef]
  36. Wang, F.; Wang, Z.; Yang, H.; Di, D.; Zhao, Y.; Liang, Q. Utilizing GRACE-Based Groundwater Drought Index for Drought Characterization and Teleconnection Factors Analysis in the North China Plain. J. Hydrol. 2020, 585, 124849. [Google Scholar] [CrossRef]
  37. Duan, Q.; Sorooshian, S.; Gupta, V. Effective and Efficient Global Optimization for Conceptual Rainfall-Runoff Models. Water Resour. Res. 1992, 28, 1015–1031. [Google Scholar] [CrossRef]
  38. Long, D.; Yang, Y.; Wada, Y.; Hong, Y.; Liang, W.; Chen, Y.; Yong, B.; Hou, A.; Wei, J.; Chen, L. Deriving Scaling Factors Using a Global Hydrological Model to Restore GRACE Total Water Storage Changes for China’s Yangtze River Basin. Remote Sens. Environ. 2015, 168, 177–193. [Google Scholar] [CrossRef]
  39. Peng, S.; Ding, Y.; Liu, W.; Li, Z. 1 Km Monthly Temperature and Precipitation Dataset for China from 1901 to 2017. Earth Syst. Sci. Data 2019, 11, 1931–1946. [Google Scholar] [CrossRef] [Green Version]
  40. Zhao, Z.; Yu, F. Atlas of Groundwater Sustainable Utilization in North China Plain; China Cartographic Publishing House: Beijing, China, 2009. (In Chinese) [Google Scholar]
  41. Zhang, C.; Duan, Q.; Yeh, P.J.-F.; Pan, Y.; Gong, H.; Gong, W.; Di, Z.; Lei, X.; Liao, W.; Huang, Z.; et al. The Effectiveness of the South-to-North Water Diversion Middle Route Project on Water Delivery and Groundwater Recovery in North China Plain. Water Resour. Res. 2020, 56, e2019WR026759. [Google Scholar] [CrossRef]
  42. Liu, R.; Zhong, B.; Li, X.; Zheng, K.; Liang, H.; Cao, J.; Yan, X.; Lyu, H. Analysis of Groundwater Changes (2003–2020) in the North China Plain Using Geodetic Measurements. J. Hydrol. Reg. Stud. 2022, 41, 101085. [Google Scholar] [CrossRef]
  43. Xu, Y.; Gong, H.; Chen, B.; Zhang, Q.; Li, Z. Long-Term and Seasonal Variation in Groundwater Storage in the North China Plain Based on GRACE. Int. J. Appl. Earth Obs. Geoinf. 2021, 104, 102560. [Google Scholar] [CrossRef]
  44. Zhao, Q.; Zhang, B.; Yao, Y.; Wu, W.; Meng, G.; Chen, Q. Geodetic and Hydrological Measurements Reveal the Recent Acceleration of Groundwater Depletion in North China Plain. J. Hydrol. 2019, 575, 1065–1072. [Google Scholar] [CrossRef]
  45. Jeong, S.-J.; Ho, C.-H.; Piao, S.; Kim, J.; Ciais, P.; Lee, Y.-B.; Jhun, J.-G.; Park, S.K. Effects of Double Cropping on Summer Climate of the North China Plain and Neighbouring Regions. Nat. Clim. Chang. 2014, 4, 615–619. [Google Scholar] [CrossRef] [Green Version]
  46. Zhang, X.; Ren, L.; Feng, W. Comparison of the Shallow Groundwater Storage Change Estimated by a Distributed Hydrological Model and GRACE Satellite Gravimetry in a Well-Irrigated Plain of the Haihe River Basin, China. J. Hydrol. 2022, 610, 127799. [Google Scholar] [CrossRef]
  47. Xiong, J.; Yin, J.; Guo, S.; Yin, W.; Rao, W.; Chao, N.; Abhishek. Using GRACE to Detect Groundwater Variation in North China Plain after South–North Water Diversion. Groundwater 2022, 61, 402–420. [Google Scholar] [CrossRef]
  48. Zhou, H.; Dai, M.; Wei, M.; Luo, Z. Quantitative Assessment of Shallow Groundwater Sustainability in North China Plain. Remote Sens. 2023, 15, 474. [Google Scholar] [CrossRef]
  49. Leng, G.; Tang, Q.; Huang, M.; Leung, L.R. A Comparative Analysis of the Impacts of Climate Change and Irrigation on Land Surface and Subsurface Hydrology in the North China Plain. Reg. Environ. Chang. 2014, 15, 251–263. [Google Scholar] [CrossRef]
  50. Zhao, A.; Xiang, K.; Zhang, A.; Zhang, X. Spatial-Temporal Evolution of Meteorological and Groundwater Droughts and Their Relationship in the North China Plain. J. Hydrol. 2022, 610, 127903. [Google Scholar] [CrossRef]
  51. Liu, J.; Jiang, L.; Zhang, X.; Druce, D.; Kittel, C.M.M.; Tøttrup, C.; Bauer-Gottwein, P. Impacts of Water Resources Management on Land Water Storage in the North China Plain: Insights from Multi-Mission Earth Observations. J. Hydrol. 2021, 603, 126933. [Google Scholar] [CrossRef]
  52. Gao, Y.; Yu, M. Assessment of the Economic Impact of South-to-North Water Diversion Project on Industrial Sectors in Beijing. J. Econ. Struct. 2018, 7, 4. [Google Scholar] [CrossRef] [Green Version]
  53. Long, D.; Yang, W.; Scanlon, B.R.; Zhao, J.; Liu, D.; Burek, P.; Pan, Y.; You, L.; Wada, Y. South-to-North Water Diversion Stabilizing Beijing’s Groundwater Levels. Nat. Commun. 2020, 11, 3665. [Google Scholar] [CrossRef]
Remotesensing 15 03264 g001 550
Figure 1. Study area of the NCP.
Figure 1. Study area of the NCP.
Remotesensing 15 03264 g001
Remotesensing 15 03264 g002 550
Figure 2. Comparison of different precipitation and actual evapotranspiration data sources in the NCP from 2003 to 2019: (a) monthly precipitation; (b) actual evapotranspiration.
Figure 2. Comparison of different precipitation and actual evapotranspiration data sources in the NCP from 2003 to 2019: (a) monthly precipitation; (b) actual evapotranspiration.
Remotesensing 15 03264 g002
Remotesensing 15 03264 g003 550
Figure 3. Distribution of parameter zones and Dirichlet boundaries (a) and geological cross section of the NCP (b) in the study area.
Figure 3. Distribution of parameter zones and Dirichlet boundaries (a) and geological cross section of the NCP (b) in the study area.
Remotesensing 15 03264 g003
Remotesensing 15 03264 g004 550
Figure 4. Simulated GWSA from January 2003 to December 2019.
Figure 4. Simulated GWSA from January 2003 to December 2019.
Remotesensing 15 03264 g004
Remotesensing 15 03264 g005 550
Figure 5. Comparison of groundwater storage changes before and after downscaling in December 2019. (a,c) The 1° and 0.05° mesh grid blocks, respectively; (b,d) the spatial distribution of GWSA before and after downscaling.
Figure 5. Comparison of groundwater storage changes before and after downscaling in December 2019. (a,c) The 1° and 0.05° mesh grid blocks, respectively; (b,d) the spatial distribution of GWSA before and after downscaling.
Remotesensing 15 03264 g005
Remotesensing 15 03264 g006 550
Figure 6. Time-series trend of downscaling results from 2003 to 2020. The red bands represent the changes in GWSA within all 0.05° grids under the 1° grid, and the blue lines represent the GWSA retrieved by GRACE and GLDAS in the 1° grid.
Figure 6. Time-series trend of downscaling results from 2003 to 2020. The red bands represent the changes in GWSA within all 0.05° grids under the 1° grid, and the blue lines represent the GWSA retrieved by GRACE and GLDAS in the 1° grid.
Remotesensing 15 03264 g006
Remotesensing 15 03264 g007 550
Figure 7. Comparison of time-series groundwater storage changes from the model and field observations.
Figure 7. Comparison of time-series groundwater storage changes from the model and field observations.
Remotesensing 15 03264 g007
Remotesensing 15 03264 g008 550
Figure 8. Comparison of TWSA of (a) Li data, (b) CLSM data, and (c) Mo data before and after standardization with the GRACE-retrieved data.
Figure 8. Comparison of TWSA of (a) Li data, (b) CLSM data, and (c) Mo data before and after standardization with the GRACE-retrieved data.
Remotesensing 15 03264 g008
Remotesensing 15 03264 g009 550
Figure 9. Change in GWSA from 2003 to 2020 (the GWSA data from July 2017 to May 2018 were calculated using Li data, Mo data, and GLDAS CLSM model data).
Figure 9. Change in GWSA from 2003 to 2020 (the GWSA data from July 2017 to May 2018 were calculated using Li data, Mo data, and GLDAS CLSM model data).
Remotesensing 15 03264 g009
Remotesensing 15 03264 g010 550
Figure 10. Change in the GWSA in the NCP and administrative regions (a) and change in GWSA and precipitation in the NCP (b) from 2003 to 2020.
Figure 10. Change in the GWSA in the NCP and administrative regions (a) and change in GWSA and precipitation in the NCP (b) from 2003 to 2020.
Remotesensing 15 03264 g010
Remotesensing 15 03264 g011 550
Figure 11. Spatial distribution map of the GWSA change rate in the NCP in four subperiods: (a) 2003–2008, (b) 2009–2014, (c) 2015–2017, and (d) 2018–2020.
Figure 11. Spatial distribution map of the GWSA change rate in the NCP in four subperiods: (a) 2003–2008, (b) 2009–2014, (c) 2015–2017, and (d) 2018–2020.
Remotesensing 15 03264 g011
Remotesensing 15 03264 g012 550
Figure 12. Spatial distribution map of the GGDI change rate: (a) 2003.1 to 2008.12, (b) 2009.1 to 2014.12, (c) 2015.1 to 2017.12, and (d) 2018.1 to 2020.12.
Figure 12. Spatial distribution map of the GGDI change rate: (a) 2003.1 to 2008.12, (b) 2009.1 to 2014.12, (c) 2015.1 to 2017.12, and (d) 2018.1 to 2020.12.
Remotesensing 15 03264 g012
Table
Table 1. The sorting of the GGDI.
Table 1. The sorting of the GGDI.
GradeClassificationGGDI
INo drought−0.5 < GGDI
IIMild drought−1.0 < GGDI ≤ −0.5
IIIModerate drought−1.5 < GGDI ≤ −1.0
IVSevere drought−2.0 < GGDI ≤ −1.5
VExtreme droughtGGDI ≤ −2.0
Table
Table 2. Data sources used in the study.
Table 2. Data sources used in the study.
Data CategoryData SourceSpatial ResolutionTime ScaleTime Span
TWSGRACE0.5°Monthly2003–2020
Li et al. [33]0.5°MonthlyJanuary 2003–June 2020
Mo et al. [34]1°Monthly2003–2020
GLDAS CLSM1°Monthly2003–2020
SMGLDAS V2.11°Monthly2003–2020
SWEGLDAS V2.11°Monthly2003–2020
PrecipitationTRMM 3B430.25°Monthly2003–2019
ERA50.25°Monthly2003–2019
PENG0.05°Monthly2003–2020
AETMOD160.05°Monthly2003–2020
ERA50.25°Monthly2003–2019
GLEAM v3.5a0.25°Monthly2003–2019
GWLIn situ observation−Monthly,
Daily		2005–2014
2018–2019
Table
Table 3. The CC, RMSE and NSE for cells in the NCP during calibration and validation periods.
Table 3. The CC, RMSE and NSE for cells in the NCP during calibration and validation periods.
Cell IDCalibration PeriodValidation PeriodCell IDCalibration PeriodValidation Period
CCRMSE (cm EWH)NSECCRMSE (cm EWH)NSECCRMSE (cm EWH)NSECCRMSE (cm EWH)NSE
G140.954.480.650.856.200.10G410.904.040.780.434.480.87
G220.953.970.760.644.690.60G420.884.850.760.304.700.85
G230.896.370.640.477.230.42G430.825.400.640.408.510.44
G240.886.090.640.207.270.42G480.965.880.820.518.160.70
G250.758.780.230.0413.78−1.75G490.874.660.700.526.370.81
G300.933.280.840.722.760.89G500.903.750.810.414.590.89
G310.903.400.790.493.570.82G560.974.600.890.809.350.69
G320.905.900.620.3111.37−0.10G570.934.990.820.550.120.65
G330.854.240.550.328.190.38G580.854.350.700.734.940.87
G390.964.250.850.428.610.52G650.942.730.900.864.390.91
G400.903.380.730.643.630.89G660.903.220.840.575.070.88
Table
Table 4. Correlation coefficients between downscaled GWSA and GWLA of twelve wells.
Table 4. Correlation coefficients between downscaled GWSA and GWLA of twelve wells.
Well NumberCCWell NumberCC
W70.51W10.52
W80.65W20.70
W90.62W30.58
W100.72W40.69
W110.67W50.60
W120.66W60.58
Table
Table 5. Comparison of the results for different time periods and research regions.
Table 5. Comparison of the results for different time periods and research regions.
Research ScholarsResearch AreaResearch PeriodChanges in GWSA in Previous ResearchChanges in GWSA in This Study
Feng et al., 2013 [5]NCP (370,000 km2)2003–2010−2.2 ± 0.3 cm/yr−0.91 cm/yr
Feng et al., 2018 [15]NCP (320,000 km2)2003–2014−7.3 ± 1.1 km3/yr−1.40 cm/yr
Liu et al., 2022 [42]NCP2003–2014−1.66 ± 0.17 cm/yr−1.40 cm/yr
2005–2016−2.21 ± 0.15 cm/yr−1.91 cm/yr
2015–2020−2.76 ± 0.55 cm/yr−2.26 cm/yr
2003–2020−2.18 ± 0.11 cm/yr−1.89 cm/yr
Gong et al., 2018 [14]NCP2003–2015−17.7 ± 1.1 mm/yr−15.2 mm/yr
HB2005–2013−14.7 ± 1.1 mm/yr−15.9 mm/yr
TJ2005–2013−20.2 ± 0.2 cm/yr−18.6 mm/yr
Xu et al., 2021 [43]NCP2003–2017−19.96 ± 3.6 cm/yr−17.15 mm/yr
Zhao et al., 2019 [44]NCP2004–mid–2016−1.7 ± 0.1 cm/yr−1.76 cm/yr
mid–2013−mid–2016−3.8 ± 0.1 cm/yr−2.47 cm/yr
Zheng et al., 2020 [30]NCP2003–2016−17.2 ± 0.8 mm/yr−16.1 mm/yr
Table
Table 6. Four-stage change rate of the downscaled GWSA in the NCP and five administrative regions.
Table 6. Four-stage change rate of the downscaled GWSA in the NCP and five administrative regions.
RegionsSlopes of GWSA in Different Periods (cm/Month)
January 2003–December 2008January 2009–December 2014January 2015–December 2017January 2018–December 2020
NCP−0.039−0.150−0.216−0.130
BJ−0.028−0.087−0.2340.095
TJ−0.068−0.141−0.3130.048
HB−0.043−0.131−0.199−0.098
HN−0.007−0.224−0.162−0.363
SD−0.040−0.175−0.245−0.154
Table
Table 7. Four-stage monthly change rate of the GGDI in the NCP and five administrative regions.
Table 7. Four-stage monthly change rate of the GGDI in the NCP and five administrative regions.
RegionsSlopes of GGDI in Different Periods
January 2003–December 2008January 2009–December 2014January 2015–December 2017January 2018–December 2020
NCP−0.0031−0.0147−0.0206−0.0117
BJ−0.0038−0.0122−0.03010.0114
TJ−0.0054−0.0117−0.02260.003
HB−0.004−0.0137−0.0191−0.0113
HN0.001−0.0192−0.0143−0.0274
SD−0.0026−0.015−0.0197−0.0113
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

Zhang, J.; Hu, L.; Sun, J.; Wang, D. Reconstructing Groundwater Storage Changes in the North China Plain Using a Numerical Model and GRACE Data. Remote Sens. 2023, 15, 3264. https://doi.org/10.3390/rs15133264

AMA Style

Zhang J, Hu L, Sun J, Wang D. Reconstructing Groundwater Storage Changes in the North China Plain Using a Numerical Model and GRACE Data. Remote Sensing. 2023; 15(13):3264. https://doi.org/10.3390/rs15133264

Chicago/Turabian Style

Zhang, Junchao, Litang Hu, Jianchong Sun, and Dao Wang. 2023. "Reconstructing Groundwater Storage Changes in the North China Plain Using a Numerical Model and GRACE Data" Remote Sensing 15, no. 13: 3264. https://doi.org/10.3390/rs15133264

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