Reconstructing Long-Term, High-Resolution Groundwater Storage Changes in the Songhua River Basin Using Supplemented GRACE and GRACE-FO Data

Abstract

The Gravity Recovery and Climate Experiment (GRACE) enables large-scale monitoring of terrestrial water storage changes, significantly contributing to hydrology and related fields. However, the coarse resolution of groundwater storage anomaly (GWSA) data limits local-scale research utilizing GRACE and GRACE-FO missions. In this study, we develop a regional downscaling model based on the linear regression relationship between GWSA and environmental variables, reducing the grid resolution of GWSA obtained from GRACE from approximately 25 km to 1 km. First, we estimate the missing values of monthly continuous terrestrial water storage anomaly (TWSA) for the period from 2003 to 2020 using interpolated multi-channel singular spectrum analysis (IMSSA). Next, we apply the water balance equation to separate GWSA from TWSA, which is provided jointly by the Global Land Data Assimilation System (GLDAS) and the distributed ecohydrological model ESSI-3. We then employ a partial least squares regression (PLSR) model to identify the most significant environmental variables related to GWSA. Precipitation (Prec), normalized difference vegetation index (NDVI), and actual evapotranspiration (AET), with variable importance in projection (VIP) values greater than 1.0, are recognized as effective variables for reconstructing long-term, high-resolution groundwater storage changes. Finally, we downscale and reconstruct the long-term (2003–2020), high-resolution (1 km × 1 km) monthly GWSA in the Songhua River Basin using fused and supplemented GRACE/GRACE-FO data, employing either geographically weighted regression (GWR) or random forest (RF) models. The results demonstrate superior performance of the GWR model (CC = 0.995, NSE = 0.989, RMSE = 2.505 mm) compared to the RF model in downscaling. The downscaled GWSA in the Songhua River Basin not only achieves high spatial resolution but also exhibits improved accuracy when compared to in situ groundwater observation records. This research enhances understanding of spatiotemporal variations in regional groundwater due to local agricultural and industrial water use, providing a scientific basis for regional water resource management.

1. Introduction

Water is fundamental to human life, economic activities, and ecosystem sustainability, with groundwater serving as a critical resource for drinking water, irrigation, and industry [1,2,3,4,5,6]. However, groundwater depletion driven by climate change and over-extraction poses significant challenges, while short-term groundwater storage fluctuations are often linked to extreme weather events such as droughts and floods [7,8,9,10]. Understanding the spatiotemporal dynamics of groundwater storage anomalies (GWSA) is essential for assessing climate impacts and ensuring sustainable water resource management [11,12].

Traditional groundwater monitoring relies on extensive well networks, but these are often costly, time-intensive, and geographically constrained, resulting in sparse and uneven data coverage [13,14,15,16]. The introduction of GRACE satellite data in 2002 revolutionized large-scale GWSA research by providing accurate terrestrial water storage anomaly (TWSA) measurements, particularly benefiting data-scarce regions [17,18,19,20,21]. GRACE data have been widely applied to investigate groundwater depletion globally, including in the Mississippi River Basin, India, and China’s Yangtze River Basin [11,22,23].

Beyond groundwater studies, GRACE has been utilized for diverse applications such as extreme hydrometeorological events, river water storage fluctuations, and glacier melting [24,25,26,27,28,29,30,31,32,33,34,35]. Despite its utility, GWSA cannot be directly retrieved from GRACE data and requires integration with complementary datasets [2,36]. For example, Feng et al. (2013) derived GWSA in North China by subtracting GLDAS soil moisture data from GRACE TWSA [37]. Similarly, Amiri et al. (2023) analyzed groundwater trends in Yazd Province, Iran, using GRACE and GLDAS components [38]. Zhang et al. (2022) combined GLDAS and independent component analysis to extract GWSA for the Haihe River Basin, achieving strong agreement with SWAT hydrological model simulations [36].

Despite significant advancements, the application of GRACE data in hydrology and water resource management faces limitations, particularly data gaps and coarse spatial resolution [39,40]. Gaps in GRACE data, notably during the 2017–2018 transition to GRACE-FO, stem from satellite sampling constraints and cumulative observations, leading to incomplete time series [41,42]. To address these gaps, researchers have employed various reconstruction techniques, including machine learning and statistical methods [42,43,44]. For instance, Sun et al. (2019) developed a deep convolutional neural network (CNN) to reconcile discrepancies between GRACE-derived and NOAH-simulated TWSA, bridging the mission transition gap [45]. Similarly, multi-channel singular spectrum analysis (MSSA) has proven effective, with Gauer et al. (2023) demonstrating its ability to reconstruct continuous GRACE data with minimal noise and improved temporal coverage [46,47,48,49].

However, coarse spatial resolution of GRACE limits its applicability for localized studies. Current products, such as those from the University of Texas Center for Space Research (CSR), offer resolutions of 0.5° to 0.25°, derived from a 3° × 3° grid, which remains inadequate for small-scale aquifer or basin analysis [50,51,52,53]. Enhancing spatial resolution of GRACE is essential for accurately capturing TWSA and GWSA variations in finer-scale regions [54,55].

Two primary approaches exist for downscaling GRACE data, i.e., dynamic and statistical methods [56,57]. Dynamic approaches integrate GRACE observations with hydrological models, leveraging physical mechanisms, but involve complex calculations sensitive to boundary conditions [58,59,60,61]. In contrast, statistical methods establish predictive relationships between coarse-resolution GRACE data and auxiliary variables, enabling finer-scale predictions [62]. These methods, while dependent on high-quality predictor variables, offer simpler and effective alternatives for improving the spatial resolution of GRACE and have been widely adopted for localized applications [50,63,64].

Yin et al. (2018) successfully downscaled GRACE-derived GWSA from 110 km to 2 km in the North China Plain by leveraging the strong correlation between evapotranspiration (ET) data and GWSA, capturing sub-grid heterogeneity [65]. Other studies have applied multiple regression, artificial neural networks, and extreme gradient boosting to reduce TWSA resolution in Michigan’s Lower Peninsula to 0.125° (120 km²) using variables such as precipitation, NDVI, snow cover, streamflow, water levels, surface temperature, soil moisture, air temperature, and ET [66]. In the Haihe River Basin, multivariate linear regression (MLR), random forest (RF), and NoahV2.1 models were used to downscale GRACE-based GWSA from 1° to 0.25°, with the RF model showing the best performance [67]. Zhang et al. (2021) applied random forest and extreme gradient boosting to downscale GRACE and GLDAS data from 1° to 0.25° to 1 km, respectively, for GWS analysis in China (2004–2016). However, the annual temporal resolution of the GWSA data limits its ability to capture monthly variations, affecting its application in hydrological research, particularly for flood and drought events [55].

The coarse resolution of GRACE data and the missing values between GRACE and GRACE-FO remain significant constraints in hydrological studies. To address this, this study generates a continuous GRACE TWSA dataset using the interpolated multi-channel singular spectrum analysis (IMSSA) method, integrating the partial least squares regression (PLSR) model, GLDAS, and the ESSI-3 ecohydrological model. By constructing a downscaling regression model, this research aims to effectively downscale GRACE-derived GWSA data. Key contributions of this paper include: (a) improved temporal consistency of GRACE TWSA through the IMSSA method; (b) enhanced quality of ESSI_GWSA data over GLDAS_GWSA, which relies on GLDAS water storage components and continuous TWSA; (c) use of the PLSR model to identify climate variables influencing water storage changes, aiding the downscaling model; and (d) development and application of geographically weighted regression (GWR) and RF downscaling models in the Songhua River Basin, producing high-quality monthly GWSA data at a 1 km resolution for 2003–2020.

2. Study Area

The Songhua River Basin, located in northeastern China (119°52′–132°31′E, 41°42′–51°38′N), spans over 550,000 km² and is one of China’s seven major river basins. As a key tributary of the Heilongjiang River, it is divided into three sub-basins: the Nenjiang River Basin (west), the Upper Songhua River Basin (south), and the Lower Songhua River Basin (northeast) [68] (Figure 1). The basin experiences a temperate monsoon climate, with cold, dry winters (temperatures often below −30 °C) and warm, humid summers influenced by maritime air currents, reaching up to 30 °C. Annual precipitation ranges from 400 to 800 mm, with the majority falling during the summer months [69,70].

Geologically, the basin is predominantly composed of Quaternary sediments, including alluvial and lacustrine deposits that form the main aquifers. These sediments, composed of sand, gravel, and clay, support both confined and unconfined aquifers, which provide critical groundwater for agricultural, industrial, and domestic use. The mountainous areas contain fractured metamorphic and igneous rock aquifers, but their contribution to regional groundwater is minimal compared to the sedimentary aquifers.

The topography of the basin is diverse, including mountains, hills, basins, and plains. The Songnen Plain, located centrally and in the southwest, is a major agricultural hub, producing commercial crops such as soybeans, corn, and rice. Groundwater irrigation is vital for sustaining agricultural production in this area [2,71].

3. Materials and Methodology

3.1. Input Datasets

This study utilized a range of datasets, including satellite observations, reanalysis products, model outputs, and ground-based data, as detailed below.

3.1.1. GRACE/GRACE-FO Data

The GRACE mission, launched in March 2002 by NASA and the German Aerospace Center (DLR), operated until June 2017, monitoring global water resources, glacier melt, sea-level changes, and Earth’s mass distribution through gravitational field variations [72,73]. Its successor, GRACE-FO, was launched in May 2018 after an 11-month delay [74].

GRACE and GRACE-FO collected high-precision monthly TWS anomaly data, processed by institutions like NASA’s JPL, GFZ, and CSR. These data, widely applied across various disciplines, are available in spherical harmonic (SHC) and Mascon solutions, with grid spatial resolutions ranging from 0.25° to 1° [27,57]. This study uses CSR RL06 Mascon data (0.25° grid resolution) to analyze GWS changes in the Songhua River Basin from 2003 to 2020. The research baseline (2005–2010) aligns with groundwater well data availability, and TWSA-derived GWSA anomalies are calculated relative to this baseline. The CSR data are publicly accessible at https://www2.csr.utexas.edu/grace/RL06_mascons.html (accessed on 10 August 2024).

3.1.2. GLDAS Model Data

GLDAS, developed by NASA GSFC and NOAA, integrates satellite and ground-based observations using advanced land surface modeling and data assimilation techniques [75]. This study uses monthly soil moisture storage (SMS), canopy water storage (CWS), and snow water equivalent (SWE) from the GLDAS Noah-v2.1 model (0.25° resolution) for 2003–2020, with data processed based on the research baseline period. It is worth noting that the GLDAS Noah model does not include the GWS component, and the daily scale product of the GLDAS CLSM includes the GWS component.

To validate the GRACE data filled using the interpolated multi-channel singular spectrum analysis method (IMSSA), the community land surface model (CLSM) from GLDAS was processed to obtain TWSA data that aligns with the new research baseline. It is important to note that CLSM TWSA does not include surface water storage components, such as lakes, reservoirs, and rivers [17]. Previous research indicates that changes in surface water within the study area can be considered negligible compared to variations in soil moisture and GWS [2,76]. Therefore, this study focuses on verifying and analyzing the trends in TWSA changes based on the filled GRACE data.

3.1.3. ESSI-3 Model Data

The ESSI-3 model, developed by Zhang Wanchang, provides monthly averages of soil moisture storage (SMS), canopy water storage (CWS), and snow water equivalent (SWE) from 2003 to 2020. As a hydrological framework with independent intellectual property rights, ESSI-3 assesses the impacts of climate change and surface changes on hydrological processes [77,78]. Its robustness in analyzing watershed hydrological components has been widely demonstrated [79,80,81,82,83].

3.1.4. Environmental Variables

Environmental variables used for downscaling GRACE data include precipitation (prec), NDVI, actual evapotranspiration (AET), land surface temperature (LST), soil moisture (SM), and air temperature (temp), selected for their significant contributions in previous studies [14,17,84].

High-resolution precipitation data (1 km) were validated against meteorological records in the Songhua River Basin [70]. NDVI data were obtained from MODIS Terra (MOD13A3, 1 km) via NASA’s LP DAAC. AET data were derived from the SSEBop model using MODIS data from the USGS FEWS NET portal (1 km) [85].

Monthly LST averages were calculated from MODIS data (MOD11A2, 1 km). The SM dataset, developed using machine learning and validated against ground observations, combines data from the China Meteorological Administration and other sources (1 km) [86,87]. Air temperature data, downscaled using the Delta method from CRU and WorldClim datasets, were verified against independent meteorological stations [88,89,90,91].

3.1.5. Ground-Based Measurements

The groundwater level borehole measurement data were obtained from the China Geological Environment Monitoring Groundwater Level Yearbook for the years 2005–2011, 2013, 2015, and 2016, with observation intervals ranging from 5 to 20 days. This study compiled 1006 records over a decade within the study area. Groundwater level anomaly (GLA) was calculated by subtracting the long-term baseline average (2005–2010) for each observation. GLA was then converted to groundwater storage anomaly (GWSA) using the specific yield (Sy) parameter from the PCR-GLOBWB model.

3.2. Methodology

This study follows a structured approach. First, the IMSSA method was applied to fill missing values in the GRACE and GRACE-FO datasets prior to estimating GWSA. Next, using the water balance equation, GWSA was isolated from the TWSA derived from GRACE data, incorporating water storage components from the GLDAS and ESSI-3 models, such as SMS, SWE, and CWS.

Subsequently, PLSR was employed to assess the contributions of selected environmental variables to the GRACE-based GWSA, with VIP scores guiding the final selection. Using the chosen variables—Prec, NDVI, and AET—a linear regression model was developed with GWR and random forest (RF) methods. This model linked GWSA to auxiliary variables at a coarse scale (0.25°).

The derived correlation was then used to predict GWSA at a finer spatial resolution of 1 km, based on high-resolution environmental variables matching this scale. The research workflow is depicted in Figure 2.

3.2.1. Partial Least Squares Regression

In the initial phase of this study, we identified environmental variables previously shown to correlate significantly with GWSA. The selected variables included Prec, NDVI, AET, LST, SM, and temp [14,17,92].

We employed the partial least squares regression (PLSR) method to assess the contributions of these environmental variables to GWSA in the study area. PLSR is a supervised learning algorithm that integrates the strengths of principal component analysis (PCA) and multivariate linear regression (MLR). It is particularly well suited for addressing high-dimensional, multicollinear, and multi-response variable regression challenges [93].

PLSR effectively incorporates information from both the target variable and the explanatory variables, allowing for an efficient analysis of the relationships between a single target variable and multiple predictors [94]. The VIP value is utilized within PLSR to evaluate the significance of each explanatory variable in relation to the target variable’s predictive capacity. A higher VIP value indicates a greater contribution of that variable to the prediction of GWSA [95]. Typically, explanatory variables with VIP values exceeding 1.0 are recognized as of greater importance to the model.

The VIP value is calculated as the weighted sum of squares of each explanatory variable across all components, as described below:

(1)$\mathrm{VI}{\mathrm{P}}_{\mathrm{j}}=\sqrt{\mathrm{p}{\displaystyle \frac{\sum _{\mathrm{a}=1}^{\mathrm{A}} \mathrm{Sa}(\mathrm{Y},{\mathrm{t}}_{\mathrm{a}})\text{×}{\mathrm{W}}_{\mathrm{ja}}^2}{\sum _{\mathrm{a}=1}^{\mathrm{A}} \mathrm{Sa}(\mathrm{Y},{\mathrm{t}}_{\mathrm{a}})}}}$

where $\mathrm{p}$ is the number of explanatory variables; $\mathrm{A}$ is the number of selected components; $\mathrm{Sa}(\mathrm{Y},{\mathrm{t}}_{\mathrm{a}})$ is the contribution of component ta to the target variable Y; and ${\mathrm{W}}_{\mathrm{ja}}$ is the weight of the explanatory variable ${\mathrm{X}}_{\mathrm{j}}$ in component ${\mathrm{t}}_{\mathrm{a}}$.

3.2.2. Interpolated Multi-Channel Singular Spectrum Analysis

Multi-channel singular spectrum analysis (MSSA) is a time series analysis method that separates signals from noise or different sources by performing singular value decomposition on multivariate time series and spatial and temporal correlations associated with channels or time series [96]. It is commonly used for denoising, trend extraction, periodic analysis, and missing value interpolation of time series data [46,47,49]. In this study, for the Songhua River Basin, before performing missing value prediction, we first perform linear interpolation on a small part of the missing value data to ensure the integrity of the time series data. Then, the interpolated MSSA method is used to predict the missing GRACE TWS data, and finally a continuously reconstructed TWS dataset is obtained. The method is mainly divided into the following steps:

• (1). Constructing the Embedding Matrix

Select a window length L (L ≪ N). For each time series, construct its trajectory matrix:

(2)${\mathrm{X}}_{\mathrm{k}}=\left[\begin{array}{cccc}{\mathrm{x}}_{\mathrm{k},1}& {\mathrm{x}}_{\mathrm{k},2}& \text{⋯}& {\mathrm{x}}_{\mathrm{k},\mathrm{N}-\mathrm{L}+1}\\ {\mathrm{x}}_{\mathrm{k},2}& {\mathrm{x}}_{\mathrm{k},3}& \text{⋯}& {\mathrm{x}}_{\mathrm{k},\mathrm{N}-\mathrm{L}+2}\\ \text{⋮}& \text{⋮}& \text{⋱}& \text{⋮}\\ {\mathrm{x}}_{\mathrm{k},\mathrm{L}}& {\mathrm{x}}_{\mathrm{k},\mathrm{L}+1}& \text{⋯}& {\mathrm{x}}_{\mathrm{k},\mathrm{N}}\end{array}\right]$

The trajectory matrix is stacked vertically to generate the embedding matrix X:

(3)$\mathrm{X}=\left[\begin{array}{c}{\mathrm{X}}_1\\ {\mathrm{X}}_2\\ \text{⋮}\\ {\mathrm{X}}_{\mathrm{K}}\end{array}\right]$

where K is the number of spatial locations and N is the length of each time series.

• (2). Singular Value Decomposition (SVD)

(4)$\mathrm{X}={\sum }_{\mathrm{i}=1}^{\mathrm{r}} {\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{U}}_{\mathrm{i}}{\mathrm{V}}_{\mathrm{i}}^{\mathrm{T}}$

where r is the rank of the embedding matrix. The singular values ${\mathsf{\lambda }}_{\mathrm{i}}$ indicate the importance of each decomposed component. High values correspond to the main trends and periodic patterns in the data, and low values usually correspond to noise. The left singular vectors ${\mathrm{U}}_{\mathrm{i}}$ represent the main patterns of the data in different channels in the time series. The right singular vectors ${\mathrm{V}}_{\mathrm{i}}$ represent the time pattern, revealing important time structures and patterns in the time series, such as seasonal changes or trends.

• (3). Time Series Reconstruction

Time series reconstruction is based on the previously calculated main components (singular values and singular vectors) to reconstruct a matrix that approximates the original time series. This matrix can fill in missing values and remove noise. For details, please refer to [97,98].

3.2.3. Groundwater Storage Anomaly Estimation

The TWSA derived from GRACE represents changes in surface water and groundwater storage relative to the long-term average observed from 2005 to 2010. This encompasses various components, including soil moisture, groundwater, snow and ice storage, as well as the water volume in lakes and rivers [17]. Previous analyses and related studies indicate that the primary contributors to changes in TWSA within the Songhua River Basin are soil moisture storage anomaly (SMSA), snow water equivalent anomaly (SWEA), and canopy water storage anomaly (CWSA) [2]. Consequently, TWSA and GWSA in this study can be calculated using the following formula:

(5)$\mathrm{TWSA}=\mathrm{GWSA}+\mathrm{SMSA}+\mathrm{SWEA}+\mathrm{CWSA}$

(6)$\mathrm{GWSA}=\mathrm{TWSA}-(\mathrm{SMSA}+\mathrm{SWEA}+\mathrm{CWSA})$

where GWSA is the groundwater storage anomaly; SMSA is the soil moisture anomaly; SWEA is the snow water equivalent anomaly; and CWSA is the canopy water storage anomaly.

3.2.4. Random Forest Method

The RF algorithm is an ensemble machine learning method based on decision trees [99]. Initially proposed by Breiman et al. (2001), it has since gained widespread application in classification, regression, feature selection, and various other research domains [100,101]. RF operates by constructing multiple decision trees, each built from a random subset of the original dataset and a randomly selected subset of features. The final prediction is derived from averaging the predictions of all individual trees. The process consists of several key steps:

Data Preparation: Match high-resolution auxiliary variables with low-resolution target variables to create input data, which is then randomly divided into training and validation sets.

Feature Selection: For each training sample in the training set, select a subset of features to construct a decision tree.

Model Fitting: Fit the regression model to the training samples and generate predictions from the individual decision tree.

Model Aggregation: Combine all individual decision trees to form a random forest model, yielding the optimal prediction.

Parameter optimization is a critical step in enhancing model performance and improving prediction accuracy. Among the parameters, the number of trees significantly influences the model’s predictive results. In this study, we employed a random grid search method to identify the optimal parameter values.

In this section, we matched the selected explanatory variables with the spatial resolution of the GWSA derived from GRACE to establish a statistical relationship between GWSA and the explanatory variables at the original resolution. Utilizing these explanatory variables at fine resolution, we constructed an RF regression model to predict GWSA at a finer scale. Finally, a residual correction was performed based on the predicted results and the original GWSA data.

3.2.5. Geographically Weighted Regression Model

Geographically weighted regression (GWR) is a spatial local regression model that not only constructs a dynamic relationship between the target variable and the explanatory variable but also introduces spatial weights into the model to deal with the spatial heterogeneity of the variables [102]. The traditional regression model assumes that the regression coefficient is constant throughout the study area; that is, the relationship between the variables is spatially consistent [103]. For example, changes in precipitation and groundwater reserves have different patterns and relationships in different spatial locations [70]. Therefore, when downscaling GWSA, the GWR method also performs well by taking into account the spatial non-stationary relationship between the dependent variable and the predictor variable [17,104]. The model can be expressed as follows:

(7)${\mathrm{y}}_{\mathrm{i}}={\mathrm{\beta }}_0\left({\mathrm{\mu }}_{\mathrm{i}},{\mathrm{\upsilon }}_{\mathrm{i}}\right)+{\sum }_{\mathrm{k}=1}^{\mathrm{n}}{\mathrm{\beta }}_{\mathrm{k}}\left({\mathrm{\mu }}_{\mathrm{i}},{\mathrm{\upsilon }}_{\mathrm{i}}\right){\mathrm{x}}_{\mathrm{ik}}+{\mathrm{\epsilon }}_{\mathrm{i}}$

where ${\mathrm{y}}_{\mathrm{i}}$ represents the dependent variable; ${\mathrm{x}}_{\mathrm{ik}}$ represents the kth independent variable; $\left({\mathrm{\mu }}_{\mathrm{i}},{\mathrm{\upsilon }}_{\mathrm{i}}\right)$ represents the geographic coordinates of the i-th point; ${\mathrm{\beta }}_0\left({\mathrm{\mu }}_{\mathrm{i}},{\mathrm{\upsilon }}_{\mathrm{i}}\right)$ represents the intercept of the i-th point; ${\mathrm{\beta }}_{\mathrm{k}}\left({\mathrm{\mu }}_{\mathrm{i}},{\mathrm{\upsilon }}_{\mathrm{i}}\right)$ represents the coefficient of ${\mathrm{x}}_{\mathrm{ik}}$; and ${\mathrm{\epsilon }}_{\mathrm{i}}$ represents the residual of the i-th point.

The regression parameters in the GWR model change with the change in spatial location, and the calculation formula of the parameters is as follows:

(8)$\mathrm{\beta }\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)={\left({\mathrm{X}}^{\mathrm{T}}\mathrm{w}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)\mathrm{X}\right)}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{w}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)\mathrm{y}$

where $\mathrm{\beta }\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$ represents the local coefficient of the position $\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$; X and y are the vectors of the explanatory variables and the dependent variables respectively; and $\mathrm{w}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$ represents the weight matrix of the i-th point.

The commonly used methods for calculating the weight matrix include the Gaussian kernel function (Gaussian), exponential kernel function (Exponential), quadratic kernel function (Bi-square), cubic kernel function (Tri cube), etc. The method used in this study is the Gaussian kernel function:

(9)${\mathrm{W}}_{\mathrm{ij}}=\exp \left(-({\mathrm{d}}_{\mathrm{ij}}/\mathrm{b}{)}^2\right)$

where ${\mathrm{W}}_{\mathrm{ij}}$ is the weight of observation position j; ${\mathrm{d}}_{\mathrm{ij}}$ is the distance between observation points i and j; and b is the bandwidth size of the kernel function.

Determining the appropriate bandwidth parameter is crucial to the accuracy of geographically weighted regression estimation. Therefore, this study selected the improved AIC information criterion (AICc) method to determine the optimal modeling bandwidth.

(10)$\mathrm{AICc}=2\mathrm{nln}\left(\widehat{\mathrm{\sigma }}\right)+\mathrm{nln}\left(2\mathrm{\pi }\right)+\mathrm{n}\left[{\displaystyle \frac{\mathrm{n}+\mathrm{tr}\left(\mathrm{S}\right)}{\mathrm{n}-2-\mathrm{tr}\left(\mathrm{S}\right)}}\right]$

where n is the number of sample points; the matrix S is the projection matrix from the observed value to the fitted value; $\mathrm{tr}\left(\mathrm{S}\right)$ represents the trajectory of the hat matrix; and $\widehat{\mathsf{\sigma }}$ is the maximum likelihood estimate of the random error term.

3.2.6. Hydrological Model ESSI-3

In this study, we employed the distributed ecohydrological model ESSI-3 to calculate monthly SMS, CWS, and SWE at a spatial resolution of 1 km. Anomalies of these variables were subsequently derived based on the baseline period from 2005 to 2010. ESSI-3 is a grid-based distributed hydrological model that simulates three layers of soil aquifers, each exhibiting distinct thicknesses and significant parameter heterogeneity. The model captures the fluxes between these layers while accounting for the energy–water interaction processes at the soil–atmosphere interface and shallow groundwater.

In addition to modeling vertical hydrological fluxes and storage dynamics, the ESSI-3 model incorporates horizontal runoff processes. Runoff is classified into three components: surface runoff (including snowmelt runoff), soil flow, and subsurface runoff. These components converge into river channels and outlets via various media, including slopes, soils, and underground aquifers. For a more detailed description of the ESSI-3 model, please refer to the related studies [77,78,79,80,81,82,105,106].

To effectively drive the ESSI-3 model, we gathered essential model-driving data on meteorological conditions, soil characteristics, vegetation types, topography, and other relevant factors. Table 1 provides details regarding the type, spatial resolution, and temporal resolution of the driving data. The ESSI-3 model was executed at a spatial resolution of 1 km for the period from 2000 to 2020. The first three years were designated as the model’s warm-up period, followed by verification of runoff simulation results using available measured hydrological data from selected years.

3.2.7. Evaluation Indicators

This study employs three indicators to assess the downscaling results and the simulation outcomes of the hydrological model: Root Mean Square Error (RMSE), Nash-Sutcliffe Efficiency (NSE), and Correlation Coefficient (CC). Both CC and NSE values range from 0 to 1, with higher values indicating superior simulation and downscaling performance. Conversely, a lower RMSE value signifies greater accuracy in the predicted GWSA. The formulas for calculating these indicators are as follows:

(11)$\mathrm{CC}={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{X}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{X}}\left)\right({\mathrm{Y}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{Y}})}{\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{X}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{X}})}^2}\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{Y}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{Y}})}^2}}}$

(12)$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}\sum _{\mathrm{i}=1}^{\mathrm{n}} {\left({\mathrm{Y}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{i}}\right)}^2}$

(13)$\mathrm{NSE}=1-{\displaystyle \frac{\sum _{\mathrm{i}=1}^{\mathrm{n}} {\left({\mathrm{Y}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{i}}\right)}^2}{\sum _{\mathrm{i}=1}^{\mathrm{n}} {\left({\mathrm{X}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{X}}\right)}^2}}$

where Y represents the observed value, X represents the predicted value; $\stackrel{\mathrm{-}}{\mathrm{Y}}$ and $\stackrel{\mathrm{-}}{\mathrm{X}}$ represent the average values of Y and X, respectively; and n represents the number of datasets.

4. Results

4.1. Reconstruction of Missing GRACE/GRACE-FO Data

During the period from 2002 to 2020, a total of 33 months of data were missing from the GRACE and GRACE-FO datasets. To solve this problem, we used interpolated multi-channel singular spectrum analysis (IMSSA), firstly filling in some missing values of GRACE TWSA data with linear interpolation, then predicting the large-scale missing data between GRACE and GRACE-FO data, and finally obtaining complete GRACE TWSA data. The CLSM of GLDAS is currently the only global model that incorporates advanced data assimilation processes, including elegant estimation, to generate hydrological fluxes [20,39,74]. We utilized this model to verify the trends observed in the reconstructed data. The results indicate a high correlation coefficient (CC) of 0.88 between the reconstructed GRACE TWSA and the CLSM TWSA (Figure 3). This strong correlation demonstrates that the changing trends of the reconstructed GRACE TWSA align closely with those of the CLSM TWSA.

To further validate the accuracy of the IMSSA method, we randomly selected a subset of observations spanning 24 months and compared the predicted TWSA with the original GRACE TWSA. The results demonstrate that the IMSSA method exhibits robust performance, achieving high accuracy. As illustrated in Figure 4, the CC between the predicted and observed results was 0.981, while the RMSE was 3.715 mm. These findings confirm the feasibility of using the IMSSA method for predicting missing GRACE TWSA data. Additionally, recent studies have further substantiated the effectiveness of the IMSSA method in reconstructing TWSA data [46,47,49].

4.2. Performance of Hydrological Models and Determination of Water Storage Components

In this study, GWSA was first isolated from GRACE TWSA by utilizing the water storage components—SMS, SWE, and CWS—derived from the GLDAS and ESSI-3 models, in accordance with the water balance equation. Subsequently, a linear regression model was established to explore the relationship between environmental variables and GWSA, which was then applied in the downscaling process.

Prior to constructing the downscaling model, we compared GWSA estimates derived from different models, specifically CLDAS_GWSA and ESSI_GWSA. The results from both GWSA models were validated using groundwater well observation data collected from the study area. Before verifying the GWSA observational data, the performance of the ESSI-3 model was evaluated. To achieve this, runoff observations from the Xiaoergou and Jiamusi hydrological stations in the Songhua River Basin were selected for calibration and validation. The model performance was assessed using the NSE and CC as evaluation metrics as shown in Table 2. The observation period (2010–2020) was divided into a calibration phase (2010–2014) and a validation phase (2015–2020). The ESSI-3 model’s performance during both phases across the two hydrological stations is summarized in Figure 5 and Table 3. The NSE values for both the calibration and validation periods at the stations were approximately 0.8, indicating strong model applicability within the Songhua River Basin. These findings demonstrate the effectiveness of the ESSI-3 model in simulating hydrological processes in the region, laying a solid foundation for the next phase of research.

To address the sparse distribution of groundwater well observation points, we divided the study area into 1° × 1° grids to ensure the continuity of data across different time periods and to improve the reliability of data verification. Each grid was assigned a unique identifier, resulting in nine grids labeled A1 to A9, as shown in Figure 6. Groundwater well observation records within each grid were processed by calculating their mean values and adjusting for the baseline period. The accuracy of these processed records was then evaluated against the GWSA data generated by various models within each grid. Table 3 summarizes the verification results, highlighting that the GWSA derived from the ESSI-3 model (ESSI_GWSA) demonstrated a stronger correlation with groundwater well observations and a lower RMSE than the GWSA derived from GLDAS (GLDAS_GWSA). Based on these findings, the ESSI_GWSA was selected for further downscaling model development.

4.3. Selection of Environmental Variables

The selection of environmental variables is critical to the construction of downscaling models [70]. After reviewing numerous studies on the downscaling of GRACE data and considering the specific conditions of the Songhua River Basin, we selected Prec, LST, AET, NDVI, Temp, and SM as potential environmental variables [14,17,62,67,92]. While human activities significantly impact GWSA variations, they are often not included in GWSA downscaling models due to the challenges related to the spatial and temporal scales of human activity data and data availability [14]. As a result, human activity factors are rarely considered as predictors in such studies.

The six candidate environmental variables, along with the GWSA separated from GRACE TWSA, were used as input and target variables, respectively, in the PLSR model. The VIP values for each environmental variable, as illustrated in Figure 7, reveal that the VIP values of Prec, AET, and NDVI are all greater than 1.0. This indicates that these three variables have the strongest influence on GWSA in the Songhua River Basin. Consequently, Prec, AET, and NDVI were selected as the final environmental variables for the GWSA downscaling model, which is based on GWR and RF methods.

4.4. Comparison of Downscaling Models

We employed Prec, AET, and NDVI as auxiliary variables, with GWSA from GRACE TWSA as the target, to develop GWR and RF downscaling models at a coarse scale (0.25°). These models were applied to produce high-resolution (1 km) downscaled results.

Both GWR and RF models have distinct advantages in downscaling applications [17,92]. GWR captures spatial heterogeneity through local regression, making it suitable for areas with complex climate and topography, such as the Songhua River Basin. It better reflects spatial trends of GWSA data but requires optimization of the geographic weight matrix and has relatively low computational efficiency. Additionally, it struggles with fitting multidimensional nonlinear relationships. In contrast, the RF model handles nonlinear relationships and complex interactions more effectively, making it ideal for high-dimensional features and climate factor modeling. However, it poorly describes spatial continuity, particularly in data-sparse areas, which may introduce some deviations. To compare model performance, we evaluated both methods using three metrics: CC, NSE, and RMSE.

The downscaling results are presented in Figure 8 and Figure 9 and Table 4. Figure 8 illustrates the annual average spatial changes in GWSA in the Songhua River Basin from 2003 to 2020, both before and after downscaling. Notably, the downscaled results, GWR_GWSA and RF_GWSA, provide more detailed spatial information compared to the coarse spatial resolution of the original GWSA while preserving the spatial distribution trends of the original data. This indicates the effectiveness of the downscaling models.

In terms of data quality, the two downscaling models demonstrate differing performances. As shown in these figures, GWR_GWSA achieves higher accuracy with a CC of 0.995, an NSE of 0.989, and n RMSE of 2.505 mm, while RF_GWSA exhibits lower accuracy, with a CC of 0.950, an NSE of 0.893, and an RMSE of 7.668 mm, respectively.

Additionally, regarding spatial distribution, RF_GWSA displays a certain degree of spatial discreteness, whereas GWR_GWSA maintains the same spatial continuity as the original GWSA data. This observation suggests that the GWR downscaling model effectively incorporates the influence of surrounding spatial information. This also proves our previous analysis of the differences between the GWR and RF models. Overall, these results indicate that the GWR downscaling model outperforms RF_GWSA in constructing a downscaling model for GWSA in the Songhua River Basin.

4.5. Analysis of GWSA Downscaling Results

Among the two downscaling schemes evaluated, the GWR model demonstrated superior performance, achieving higher NSE values and lower RMSE compared to the alternative approach. To further elucidate the downscaling results of GWR_GWSA, we compared the GWSA at the original grid resolution of 0.25° with the downscaled GWSA at 1 km.

As illustrated in Figure 10 and Figure 11, the downscaled GWSA exhibits a spatial distribution that is consistent with the original GWSA derived from GRACE. Both datasets reveal higher spatial distributions in the central and southern regions of the study area, with lower values observed in the eastern and western areas. The spatial distribution maps of GWSA before and after downscaling across different months demonstrate a strong correspondence between the two datasets. This indicates that the downscaling model effectively preserves the original spatial variation information while refining the GWSA data to a resolution of 1 km.

Figure 12 and Figure 13 illustrate that, compared to the data quality prior to downscaling, the results obtained after downscaling exhibit enhanced spatial resolution of GWSA and improved data accuracy. This improvement further validates the effectiveness of the GWR model for GWSA downscaling.

Figure 14 illustrates the geographical changes in GWSA in 2005, comparing GRACE-derived values with GWR-scaled values. The high-resolution GWSA measurements along the sample lines reveal intricate spatial variations that are obscured in the coarse-resolution data. The three sample lines (L1, L2, and L3) encompass 12, 18, and 20 pixels of the original GWSA data, respectively (Figure 14c), while the downscaled GWSA data correspond to 337, 521, and 568 pixels (Figure 14d). This downscaling effectively captures local spatial variations that are challenging to detect at the original resolution. Notably, the downscaled data exhibit significant spatial variability across the sample lines, whereas the original GWSA is represented by homogeneous pixels, as shown in Figure 14c. This indicates that the downscaled results not only align with the original GWSA but also provide enhanced detail in spatial heterogeneity.

5. Discussion

5.1. Performance of the Proposed Downscaling Model and Method

The GWR downscaling method, which integrates GRACE TWSA data with the ESSI-3 model, demonstrates robust performance in reconstructing high-resolution GWSA in the Songhua River Basin. While the GLDAS is commonly employed for downscaling GRACE data, this study verifies and compares GWSA derived from both the ESSI-3 model and GLDAS against groundwater observation data. The results indicate that the ESSI-3 model not only exhibits strong hydrological simulation performance (NSE > 0.8) but also shows a higher consistency with observed data. Utilizing this model as a target variable enhances the quality of the downscaling results.

The downscaling model’s effectiveness is fundamentally based on the linear regression relationship between the target and environmental variables. Consequently, the selection of environmental variables is critical for achieving accurate downscaling outcomes. In this study, the PLSR method was employed to assess the importance of various environmental variables, identifying those with VIP values greater than 1.0, specifically, Prec, AET, and NDVI, as the final environmental variables.

Moreover, the choice of downscaling model significantly influences the results. A comparative analysis of the GWR and RF models revealed that the GWR model exhibited more stable performance than the RF model. This finding reinforces the potential for integrating the ESSI-3 model with the GWR approach for effective GRACE downscaling.

5.2. Analysis of GWSA Change Trend in the Songhua River Basin

The Songhua River Basin, located in the northeasternmost region of China, serves as a crucial freshwater source for supporting the rapid industrial and agricultural development in the area. In recent years, accelerated economic growth and climate change have exacerbated the depletion of water resources, leading to insufficient groundwater supplies in certain regions and even potential exhaustion [2,51,68,107,108]. Consequently, monitoring the long-term, high-resolution spatiotemporal changes in groundwater is vital for sustainable water management in the basin.

Figure 15 and Figure 16 illustrate the long-term changes in GWSA values and the spatial distribution of trend values in the Songhua River Basin, respectively. As depicted in Figure 15, the fluctuations in both the original monthly GWSA and the downscaled GWSA data are largely consistent throughout the study period, further demonstrating the effectiveness of the GWR model for GWSA downscaling.

From 2003 to 2020, the downscaled GWSA exhibited an overall declining trend at a rate of −0.50 ± 0.42 mm/year, while the original GWSA decreased at a rate of −0.49 ± 0.41 mm/year. The monthly GWSA fluctuations during this period can be categorized into four distinct phases: (I) January 2003 to July 2009, (II) August 2009 to May 2012, (III) June 2012 to April 2019, and (IV) May 2019 to December 2020, with each phase represented by different colors in Figure 15. Notably, the GWSA trend during period I was relatively stable, while periods II and III exhibited a clear upward trajectory, followed by a decline in period IV.

The annual GWSA trend in the Songhua River Basin can also be divided into three periods: from 2003 to 2008, the GWSA increased at a rate of 1.04 mm/year; from 2009 to 2013, it decreased at −2.39 mm/year; and from 2014 to 2020, it further declined at a rate of −0.68 mm/year. These findings align with the research results of other scholars studying GWSA in the Songhua River Basin [2].

The spatial distribution of GWSA trends post-downscaling reveals that most areas in the eastern, southern, and western regions of the Songhua River Basin exhibit a downward trend, while some northern and northeastern areas show an upward trend. This spatial variation is primarily attributed to urban development and agricultural irrigation, as indicated by land use data.

5.3. Analysis of LULC and GWSA Changes

Land use/land cover (LULC) is one of the indicators that directly reflect changes in human activities. This section discusses the response relationship between human activities and GWSA changes by analyzing the spatiotemporal changes of LULC and GWSA from 2003 to 2020.

As shown in Figure 16a and Figure 17, and Table 5, three of the six major land use types exhibited a decline during the study period, with the water body experiencing the largest decrease (−26%). The remaining three types showed an increase, with urban land seeing the largest rise (21%). Overall, the GWSA in the Songhua River Basin from 2003 to 2020 displayed a downward trend, closely linked to changes in land use types. The expansion of cultivated land, coupled with the reduction in the number of forest, grassland, and water bodies, likely increased agricultural irrigation demand, decreased soil water storage, and limited groundwater recharge, contributing to the decline in GWSA. The growth of urban and unused land further intensified groundwater extraction, reduced surface permeability, and hindered groundwater recharge, exacerbating the GWSA decline.

Spatially, the central and southern regions of the basin showed a more pronounced decrease in GWSA, coinciding with significant increases in urban and unused land and declines in forest, grassland, and water bodies. These land use changes have a marked impact on the downward trend of GWSA. The expansion of cultivated land and urbanization are key drivers of GWSA reduction, while the loss of forests, grasslands, and water bodies, alongside the growth of unused land, indirectly reduces groundwater recharge, further depleting groundwater resources. This analysis provides valuable insights into the mechanisms by which land use changes affect regional groundwater reserves and offers a scientific basis for water resource management and land use planning.

5.4. Limitations and Research Prospects

Although the GWSA downscaling study has yielded promising results in the Songhua River Basin, several uncertainties and challenges remain that may affect the quality of the final downscaling outputs. Firstly, to mitigate uncertainties associated with GRACE data, this study utilized mascon data, which is generally regarded as superior to the original spherical harmonic data. However, variations in the models and data processing techniques employed can introduce inherent uncertainties in the mascon solutions [109]. Future research could explore the fusion of various GRACE data products using machine learning and artificial intelligence methodologies to enhance the quality of TWSA data, thereby reducing uncertainties related to GRACE measurements [110,111].

Secondly, while an IMSSA method was successfully employed to address missing data, uncertainties persisted. Future studies could investigate alternative methods for deriving TWSA data, such as Seasonal Trend Decomposition using LOESS (STL) and other machine learning approaches, potentially integrating or comparing these results with continuous TWSA data derived from the IMSSA method.

Thirdly, uncertainties in water storage components derived from hydrological models, including the ESSI-3 and GLDAS models, present additional challenges. A potential solution lies in coupling outputs from different models [112], yet the integration of hydrological models with varying research backgrounds and scales remains a recognized difficulty in hydrology.

Lastly, the limited availability of groundwater observation data in the study area affects the calibration and validation of research findings. The existing groundwater well observation data in the Songhua River Basin are relatively sparse in both temporal and spatial dimensions. Consequently, when evaluating GWSA results, we can only perform a rough verification of the observation point coverage using a grid approach. Future efforts should focus on collecting more comprehensive groundwater observation data to facilitate better analysis and validation of GRACE products. Additionally, integrating high-accuracy measurement data with downscaling results could further enhance the precision of the findings.

6. Conclusions

The coarse resolution of GRACE observations presents challenges in studying the dynamic changes in water resources at local scales. While various methods, including machine learning and neural networks, have been applied to downscale GWSA data, most studies have not addressed the estimation of missing data between GRACE and GRACE-FO prior to downscaling. Furthermore, there is a scarcity of research detailing the identification of environmental explanatory variables and the effectiveness of different downscaling schemes during model construction.

This study integrates continuous GRACE TWSA data derived from the IMSSA method with outputs from the ESSI-3 model and the GLDAS model, employing the water balance equation to generate long-term continuous GWSA data spanning from 2003 to 2020. The PLSR model identifies key climate factors influencing GWSA in the Songhua River Basin, including Prec, AET, and NDVI, with VIP scores exceeding 1.0. These factors are utilized as environmental explanatory variables for the GWSA downscaling model.

Subsequently, GWR and RF machine learning algorithms are applied to the downscaling model, successfully enhancing the spatial grid resolution of GWSA data from 0.25° to 1 km. The results indicate that the GWR model outperforms the RF model, effectively identifying spatial variations in the original GRACE-derived GWSA while preserving the overall characteristics of the data. Verification against GWSA data before and after downscaling, as well as groundwater observation records, corroborates these findings. Specifically, the CC and RMSE of the downscaled results improved by 25.6% and decreased by 14.6%, respectively, compared to the original GWSA.

In summary, the downscaling scheme proposed in this study not only captures detailed spatial variation information but also enhances data continuity. This method demonstrates significant application potential, with higher-resolution GWSA data contributing positively to the understanding of spatiotemporal changes in local water resources. Furthermore, it provides valuable quantitative information for effective regional management of agricultural and industrial water resources.

Author Contributions

Conceptualization, W.Z. and C.L.; methodology, Z.Z., C.L. and W.Z.; software, Z.Z. and C.X.; validation, C.L.; data curation, C.L. and Z.Z.; writing—original draft preparation, C.L. and W.Z.; writing—review and editing, Z.Z. and W.Z.; supervision, W.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

The authors are grateful to the Space–Time Tripolar Environmental Big Data Platform (https://doi.org/10.11888/Terre.tpdc.272415) (accessed on 20 August 2024), the National Aeronautics and Space Administration (NASA), the United States Geological Survey (https://www.usgs.gov/) (accessed on 20 August 2024), and the China Meteorological Administration (CMA) (https://data.cma.cn) (accessed on 20 August 2024) for providing datasets.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Footnotes

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.

Figures and Tables

Figure 1. Songhua River Basin: (a) Digital elevation model (DEM) along with river systems and borehole in situ points; (b) geo-location map of study site in China; (c) annual mean precipitation map from 2003 to 2020; and (d); and use and land cover map.

Figure 2. Framework for constructing a downscaling model of GWSA changes in the Songhua River Basin. (a,b): Scatter plot analysis of downscaled GWSA of different models and the original GWSA; (c–f): Spatial distribution of GWSA change trends after downscaling in the Songhua River Basin in different periods.

Figure 3. Comparison of CLSM TWSA and GRACE TWSA supplemented by the IMSSA method.

Figure 4. Evaluation and comparison of TWSA predicted by IMSSA with GRACE-TWSA using 24 months of test data from 2003 to 2020.

Figure 5. Comparison of measured and simulated discharges during the calibration period (2010–2014) and validation period (2015–2020) of two hydrological stations: (a) Xiaoergou station and (b) Jiamusi station.

Figure 6. In situ grids distribution map of the Songhua River Basin: A1–A9: In situ grid number.

Figure 7. VIP values of all variables based on the PLSR model.

Figure 8. Comparison of the spatial distribution of annual average GWSA results of different downscaling schemes and original GWSA data: (a) original GWSA; (b) GWR model results; and (c) RF model results.

Figure 9. Scatter plot analysis of downscaled GWSA of different models and the original GWSA: (a) GWR model and (b) RF model.

Figure 10. Spatial distribution of multi-year monthly average GWSA at original grid resolution in the Songhua River Basin (0.25°).

Figure 11. Spatial distribution of multi-year monthly average GWSA after downscaling in the Songhua River Basin (1 km).

Figure 12. Radar chart of the comparison results between GWSA and in situ grids before and after downscaling: (a) CC and (b) RMSE. O is the GWSA data before downscaling, and D is the GWSA data after downscaling.

Figure 13. Comparative analysis of GWR_GWSA and Insitu_GWSA measurements.

Figure 14. Spatial distribution of GWSA at different grid resolutions in 2005 and the changes with geographical location under different example lines: (a,c) 0.25° and (b,d) 1 km.

Figure 15. Time series of GRACE-derived GWSA and downscaled GWSA changes in the Songhua River Basin from 2003 to 2020: (I) January 2003 to July 2009, (II) August 2009 to May 2012, (III) June 2012 to April 2019, and (IV) May 2019 to December 2020.

Figure 16. Spatial distribution of GWSA change trends after downscaling in the Songhua River Basin in different periods: (a) 2003–2020; (b) 2003–2008; (c) 2009–2013; and (d) 2014–2020.

Figure 17. Spatial distribution of LULC and GWSA in the Songhua River Basin at different periods: (a) LULC_2003; (b) LULC_2020; (c) GWSA_2003; and (d) GWSA_2020.

Table 1

Detailed information of ESSI-3 model-driving data.

Type	Variable	Data Source	Resolution and Time
Meteorology	Precipitation	ERA5_land(https://www.ecmwf.int/en/forecasts/datasets) (accessed on 15 August 2024)	0.1°, 1 h, 2000–2020
	temperature
	Wind speed
	Surface air pressure
	Surface net solarradiation
	Surface solar radiation
	Relative humidity
Soil property	Bulk density	SoilGrids(https://www.soilgrids.org) (accessed on 15 August 2024)	1 km, fixed
	Clay content massfraction
	Silt content massfraction
	Sand content massfraction
Vegetationparameter	Leaf area index (LAI)	GLOBMAP LAI(https://zenodo.org/) (accessed on 15 August 2024)	8 km, 8-day,2000–2020
	Land use/cover(LULC)	Resources and Environment Data Cloud Platform(http://www.resdc.cn) (accessed on 15 August 2024)	1 km, yearly,2001–2020
	Tree cover fraction	MODIS(https://lpdaac.usgs.gov) (accessed on 15 August 2024)	500 m, fixed, 2010
Others	DEM	SRTMDEM(https://www.gscloud.cn) (accessed on 15 August 2024)Water Yearbook	90 m, fixedmonthly, 2010–2020
	Streamflow at Jiamusi and XiaoergouStations

Table 2

Accuracy analysis of NSE and CC at two hydrological stations during calibration (2010–2014) and validation (2015–2020) periods.

Hydrological Stations	Calibration (2010–2015)		Validation (2016–2020)		Entire Period(2010–2020)
	NSE	CC	NSE	CC	NSE	CC
Xiaoergou	0.82	0.91	0.79	0.89	0.81	0.90
Jiamusi	0.85	0.93	0.84	0.92	0.84	0.93

Table 3

Comparative analysis of the accuracy of ESSI_GWSA and GLDAS_GWSA based on in situ grids.

Grids	GLDAS GWSA		ESSI GWSA
	CC	RMSE	CC	RMSE
A1	0.106	44.811	0.509	25.561
A2	0.420	62.641	0.427	57.469
A3	0.414	85.754	0.464	67.216
A4	0.502	81.472	0.590	41.714
A5	−0.112	153.509	0.467	134.324
A6	−0.191	49.167	0.664	23.417
A7	0.332	73.991	0.218	53.525
A8	−0.313	73.664	0.484	35.731
A9	0.248	66.622	−0.315	59.016

Table 4

Comparative analysis of the accuracy of GWR_GWSA and RF_GWSA based on original GWSA data.

	CC	NSE	RMSE
GWR_GWSA	0.995	0.989	2.505
RF_GWSA	0.950	0.893	7.668

Table 5

LULC change rate in the Songhua River Basin from 2003 to 2020.

	2003 (km2)	2020 (km2)	Change Rate (%)
Crop	201,245	209,316	4
Forest	218,133	214,339	−2
Grassland	66,556	59,437	−11
Water	14,341	10,545	−26
Building	13,644	16,462	21
Unutilized land	35,936	39,500	10

References

1. Arshad, A.; Zhang, Z.; Zhang, W.; Dilawar, A. Mapping Favorable Groundwater Potential Recharge Zones Using a GIS-Based Analytical Hierarchical Process and Probability Frequency Ratio Model: A Case Study from an Agro-Urban Region of Pakistan. Geosci. Front.; 2020; 11, pp. 1805-1819. [DOI: https://dx.doi.org/10.1016/j.gsf.2019.12.013]

2. Chen, H.; Zhang, W.; Nie, N.; Guo, Y. Long-Term Groundwater Storage Variations Estimated in the Songhua River Basin by Using GRACE Products, Land Surface Models, and in-Situ Observations. Sci. Total Environ.; 2019; 649, pp. 372-387. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.08.352] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30176450]

3. Famiglietti, J.S. The Global Groundwater Crisis. Nat. Clim. Chang.; 2014; 4, pp. 945-948. [DOI: https://dx.doi.org/10.1038/nclimate2425]

4. Yari, A.; Ardalan, A.; Ostadtaghizadeh, A.; Zarezadeh, Y.; Boubakran, M.S.; Bidarpoor, F.; Rahimiforoushani, A. Underlying Factors Affecting Death Due to Flood in Iran: A Qualitative Content Analysis. Int. J. Disaster Risk Reduct.; 2019; 40, 101258. [DOI: https://dx.doi.org/10.1016/j.ijdrr.2019.101258]

5. Cuthbert, M.O.; Gleeson, T.; Moosdorf, N.; Befus, K.M.; Schneider, A.; Hartmann, J.; Lehner, B. Global Patterns and Dynamics of Climate–Groundwater Interactions. Nat. Clim. Chang.; 2019; 9, pp. 137-141. [DOI: https://dx.doi.org/10.1038/s41558-018-0386-4]

6. 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, pp. 1065-1072. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2019.06.016]

7. Ostad-Ali-Askari, K.; Shayannejad, M. Quantity and Quality Modelling of Groundwater to Manage Water Resources in Isfahan-Borkhar Aquifer. Environ. Dev. Sustain.; 2021; 23, pp. 15943-15959. [DOI: https://dx.doi.org/10.1007/s10668-021-01323-1]

8. Rodell, M.; Famiglietti, J.S.; Wiese, D.N.; Reager, J.T.; Beaudoing, H.K.; Landerer, F.W.; Lo, M.-H. Emerging Trends in Global Freshwater Availability. Nature; 2018; 557, pp. 651-659. [DOI: https://dx.doi.org/10.1038/s41586-018-0123-1]

9. Khan, A.A.; Zhao, Y.; Khan, J.; Rahman, G.; Rafiq, M.; Moazzam, M.F.U. Spatial and Temporal Analysis of Rainfall and Drought Condition in Southwest Xinjiang in Northwest China, Using Various Climate Indices. Earth Syst. Environ.; 2021; 5, pp. 201-216. [DOI: https://dx.doi.org/10.1007/s41748-021-00226-5]

10. Valipour, M.; Bateni, S.M.; Jun, C. Global Surface Temperature: A New Insight. Climate; 2021; 9, 81. [DOI: https://dx.doi.org/10.3390/cli9050081]

11. Rodell, M.; Chen, J.; Kato, H.; Famiglietti, J.S.; Nigro, J.; Wilson, C.R. Estimating Groundwater Storage Changes in the Mississippi River Basin (USA) Using GRACE. Hydrogeol. J.; 2007; 15, pp. 159-166. [DOI: https://dx.doi.org/10.1007/s10040-006-0103-7]

12. Tapley, B.D.; Watkins, M.M.; Flechtner, F.; Reigber, C.; Bettadpur, S.; Rodell, M.; Sasgen, I.; Famiglietti, J.S.; Landerer, F.W.; Chambers, D.P. et al. Contributions of GRACE to Understanding Climate Change. Nat. Clim. Chang.; 2019; 9, pp. 358-369. [DOI: https://dx.doi.org/10.1038/s41558-019-0456-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31534490]

13. Famiglietti, J.S.; Lo, M.; Ho, S.L.; Bethune, J.; Anderson, K.J.; Syed, T.H.; Swenson, S.C.; de Linage, C.R.; Rodell, M. Satellites Measure Recent Rates of Groundwater Depletion in California’s Central Valley. Geophys. Res. Lett.; 2011; 38. [DOI: https://dx.doi.org/10.1029/2010GL046442]

14. Yin, W.; Zhang, G.; Han, S.-C.; Yeo, I.-Y.; Zhang, M. Improving the Resolution of GRACE-Based Water Storage Estimates Based on Machine Learning Downscaling Schemes. J. Hydrol.; 2022; 613, 128447. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2022.128447]

15. Bierkens, M.F.P.; Wada, Y. Non-Renewable Groundwater Use and Groundwater Depletion: A Review. Environ. Res. Lett.; 2019; 14, 063002. [DOI: https://dx.doi.org/10.1088/1748-9326/ab1a5f]

16. Chen, J.; Famigliett, J.S.; Scanlon, B.R.; Rodell, M. Groundwater Storage Changes: Present Status from GRACE Observations. Surv. Geophys.; 2016; 37, pp. 397-417. [DOI: https://dx.doi.org/10.1007/s10712-015-9332-4]

17. Ali, S.; Ran, J.; Luan, Y.; Khorrami, B.; Xiao, Y.; Tangdamrongsub, N. The GWR Model-Based Regional Downscaling of GRACE/GRACE-FO Derived Groundwater Storage to Investigate Local-Scale Variations in the North China Plain. Sci. Total Environ.; 2024; 908, 168239. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2023.168239] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37931810]

18. Yang, X.; Tian, S.; Feng, W.; Ran, J.; You, W.; Jiang, Z.; Gong, X. Spatio-Temporal Evaluation of Water Storage Trends from Hydrological Models over Australia Using GRACE Mascon Solutions. Remote Sens.; 2020; 12, 3578. [DOI: https://dx.doi.org/10.3390/rs12213578]

19. Ali, S.; Liu, D.; Fu, Q.; Cheema, M.J.M.; Pal, S.C.; Arshad, A.; Pham, Q.B.; Zhang, L. Constructing High-Resolution Groundwater Drought at Spatio-Temporal Scale Using GRACE Satellite Data Based on Machine Learning in the Indus Basin. J. Hydrol.; 2022; 612, 128295. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2022.128295]

20. Khorrami, B.; Ali, S.; Gündüz, O. Investigating the Local-Scale Fluctuations of Groundwater Storage by Using Downscaled GRACE/GRACE-FO JPL Mascon Product Based on Machine Learning (ML) Algorithm. Water Resour. Manag.; 2023; 37, pp. 3439-3456. [DOI: https://dx.doi.org/10.1007/s11269-023-03509-w]

21. Ran, J.; Ditmar, P.; Liu, L.; Xiao, Y.; Klees, R.; Tang, X. Analysis and Mitigation of Biases in Greenland Ice Sheet Mass Balance Trend Estimates From GRACE Mascon Products. J. Geophys. Res. Solid Earth; 2021; 126, e2020JB020880. [DOI: https://dx.doi.org/10.1029/2020JB020880]

22. Rodell, M.; Velicogna, I.; Famiglietti, J.S. Satellite-Based Estimates of Groundwater Depletion in India. Nature; 2009; 460, pp. 999-1002. [DOI: https://dx.doi.org/10.1038/nature08238] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19675570]

23. 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, pp. 177-193. [DOI: https://dx.doi.org/10.1016/j.rse.2015.07.003]

24. 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, pp. 1791-1799. [DOI: https://dx.doi.org/10.1002/2014GL062498]

25. Tangdamrongsub, N.; Han, S.-C.; Tian, S.; Müller Schmied, H.; Sutanudjaja, E.H.; Ran, J.; Feng, W. Evaluation of Groundwater Storage Variations Estimated from GRACE Data Assimilation and State-of-the-Art Land Surface Models in Australia and the North China Plain. Remote Sens.; 2018; 10, 483. [DOI: https://dx.doi.org/10.3390/rs10030483]

26. Nourani, V.; Jabbarian Paknezhad, N.; Ng, A.; Wen, Z.; Dabrowska, D.; Üzelaltınbulat, S. Application of the Machine Learning Methods for GRACE Data Based Groundwater Modeling, a Systematic Review. Groundw. Sustain. Dev.; 2024; 25, 101113. [DOI: https://dx.doi.org/10.1016/j.gsd.2024.101113]

27. Liu, D.; Mishra, A.K.; Yu, Z.; Lü, H.; Li, Y. Support Vector Machine and Data Assimilation Framework for Groundwater Level Forecasting Using GRACE Satellite Data. J. Hydrol.; 2021; 603, 126929. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.126929]

28. Lopez, T.; Al Bitar, A.; Biancamaria, S.; Güntner, A.; Jäggi, A. On the Use of Satellite Remote Sensing to Detect Floods and Droughts at Large Scales. Surv. Geophys.; 2020; 41, pp. 1461-1487. [DOI: https://dx.doi.org/10.1007/s10712-020-09618-0]

29. Yin, J.; Gentine, P.; Slater, L.; Gu, L.; Pokhrel, Y.; Hanasaki, N.; Guo, S.; Xiong, L.; Schlenker, W. Future Socio-Ecosystem Productivity Threatened by Compound Drought–Heatwave Events. Nat. Sustain.; 2023; 6, pp. 259-272. [DOI: https://dx.doi.org/10.1038/s41893-022-01024-1]

30. Frappart, F.; Papa, F.; Famiglietti, J.S.; Prigent, C.; Rossow, W.B.; Seyler, F. Interannual Variations of River Water Storage from a Multiple Satellite Approach: A Case Study for the Rio Negro River Basin. J. Geophys. Res. Atmos.; 2008; 113, [DOI: https://dx.doi.org/10.1029/2007JD009438]

31. Hasan, E.; Tarhule, A.; Kirstetter, P.-E. Twentieth and Twenty-First Century Water Storage Changes in the Nile River Basin from GRACE/GRACE-FO and Modeling. Remote Sens.; 2021; 13, 953. [DOI: https://dx.doi.org/10.3390/rs13050953]

32. Bolaños, S.; Salazar, J.F.; Betancur, T.; Werner, M. GRACE Reveals Depletion of Water Storage in Northwestern South America between ENSO Extremes. J. Hydrol.; 2021; 596, 125687. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2020.125687]

33. Cazenave, A.; Chen, J. Time-Variable Gravity from Space and Present-Day Mass Redistribution in theEarth System. Earth Planet. Sci. Lett.; 2010; 298, pp. 263-274. [DOI: https://dx.doi.org/10.1016/j.epsl.2010.07.035]

34. Chen, J.L.; Wilson, C.R.; Tapley, B.D. Contribution of Ice Sheet and Mountain Glacier Melt to Recent Sea Level Rise. Nat. Geosci.; 2013; 6, pp. 549-552. [DOI: https://dx.doi.org/10.1038/ngeo1829]

35. Liu, B.; Zou, X.; Yi, S.; Sneeuw, N.; Li, J.; Cai, J. Reconstructing GRACE-like Time Series of High Mountain Glacier Mass Anomalies. Remote Sens. Environ.; 2022; 280, 113177. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113177]

36. 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. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2022.127799]

37. 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, pp. 2110-2118. [DOI: https://dx.doi.org/10.1002/wrcr.20192]

38. Amiri, V.; Ali, S.; Sohrabi, N. Estimating the Spatio-Temporal Assessment of GRACE/GRACE-FO Derived Groundwater Storage Depletion and Validation with in-Situ Water Quality Data (Yazd Province, Central Iran). J. Hydrol.; 2023; 620, 129416. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2023.129416]

39. Khorrami, B.; Gunduz, O. An Enhanced Water Storage Deficit Index (EWSDI) for Drought Detection Using GRACE Gravity Estimates. J. Hydrol.; 2021; 603, 126812. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.126812]

40. 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. [DOI: https://dx.doi.org/10.1029/2021GL093492]

41. Ghobadi-Far, K.; Han, S.-C.; Allgeyer, S.; Tregoning, P.; Sauber, J.; Behzadpour, S.; Mayer-Gürr, T.; Sneeuw, N.; Okal, E. GRACE Gravitational Measurements of Tsunamis after the 2004, 2010, and 2011 Great Earthquakes. J. Geod.; 2020; 94, 65. [DOI: https://dx.doi.org/10.1007/s00190-020-01395-3]

42. Sun, Z.; Long, D.; Yang, W.; Li, X.; Pan, Y. Reconstruction of GRACE Data on Changes in Total Water Storage Over the Global Land Surface and 60 Basins. Water Resour. Res.; 2020; 56, e2019WR026250. [DOI: https://dx.doi.org/10.1029/2019WR026250]

43. 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. [DOI: https://dx.doi.org/10.1029/2019WR026551]

44. Mukherjee, A.; Ramachandran, P. Prediction of GWL with the Help of GRACE TWS for Unevenly Spaced Time Series Data in India: Analysis of Comparative Performances of SVR, ANN and LRM. J. Hydrol.; 2018; 558, pp. 647-658. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2018.02.005]

45. Sun, A.Y.; Scanlon, B.R.; Zhang, Z.; Walling, D.; Bhanja, S.N.; Mukherjee, A.; Zhong, Z. Combining Physically Based Modeling and Deep Learning for Fusing GRACE Satellite Data: Can We Learn From Mismatch?. Water Resour. Res.; 2019; 55, pp. 1179-1195. [DOI: https://dx.doi.org/10.1029/2018WR023333]

46. Lai, Y.; Zhang, B.; Yao, Y.; Liu, L.; Yan, X.; He, Y.; Ou, S. Reconstructing the Data Gap between GRACE and GRACE Follow-on at the Basin Scale Using Artificial Neural Network. Sci. Total Environ.; 2022; 823, 153770. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2022.153770]

47. Wang, F.; Shen, Y.; Chen, Q.; Wang, W. Bridging the Gap between GRACE and GRACE Follow-on Monthly Gravity Field Solutions Using Improved Multichannel Singular Spectrum Analysis. J. Hydrol.; 2021; 594, 125972. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.125972]

48. Yi, S.; Sneeuw, N. Filling the Data Gaps Within GRACE Missions Using Singular Spectrum Analysis. J. Geophys. Res. Solid Earth; 2021; 126, e2020JB021227. [DOI: https://dx.doi.org/10.1029/2020JB021227]

49. Gauer, L.-M.; Chanard, K.; Fleitout, L. Data-Driven Gap Filling and Spatio-Temporal Filtering of the GRACE and GRACE-FO Records. J. Geophys. Res. Solid Earth; 2023; 128, e2022JB025561. [DOI: https://dx.doi.org/10.1029/2022JB025561]

50. Chen, Z.; Zheng, W.; Yin, W.; Li, X.; Zhang, G.; Zhang, J. Improving the Spatial Resolution of GRACE-Derived Terrestrial Water Storage Changes in Small Areas Using the Machine Learning Spatial Downscaling Method. Remote Sens.; 2021; 13, 4760. [DOI: https://dx.doi.org/10.3390/rs13234760]

51. Wang, S.; Cui, G.; Li, X.; Liu, Y.; Li, X.; Tong, S.; Zhang, M. GRACE Satellite-Based Analysis of Spatiotemporal Evolution and Driving Factors of Groundwater Storage in the Black Soil Region of Northeast China. Remote Sens.; 2023; 15, 704. [DOI: https://dx.doi.org/10.3390/rs15030704]

52. Watkins, M.M.; Wiese, D.N.; Yuan, D.-N.; Boening, C.; Landerer, F.W. Improved Methods for Observing Earth’s Time Variable Mass Distribution with GRACE Using Spherical Cap Mascons. J. Geophys. Res. Solid Earth; 2015; 120, pp. 2648-2671. [DOI: https://dx.doi.org/10.1002/2014JB011547]

53. 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. [DOI: https://dx.doi.org/10.3390/rs12132151]

54. Akhtar, F.; Nawaz, R.A.; Hafeez, M.; Awan, U.K.; Borgemeister, C.; Tischbein, B. Evaluation of GRACE Derived Groundwater Storage Changes in Different Agro-Ecological Zones of the Indus Basin. J. Hydrol.; 2022; 605, 127369. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.127369]

55. Zhang, J.; Liu, K.; Wang, M. Downscaling Groundwater Storage Data in China to a 1-Km Resolution Using Machine Learning Methods. Remote Sens.; 2021; 13, 523. [DOI: https://dx.doi.org/10.3390/rs13030523]

56. Atkinson, P.M. Downscaling in Remote Sensing. Int. J. Appl. Earth Obs. Geoinf.; 2013; 22, pp. 106-114. [DOI: https://dx.doi.org/10.1016/j.jag.2012.04.012]

57. Yazdian, H.; Salmani-Dehaghi, N.; Alijanian, M. A Spatially Promoted SVM Model for GRACE Downscaling: Using Ground and Satellite-Based Datasets. J. Hydrol.; 2023; 626, 130214. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2023.130214]

58. Fowler, H.J.; Blenkinsop, S.; Tebaldi, C. Linking Climate c Hange Modelling t o Impacts Studies: Recent advances in downscaling t echniques f or hydrological modelling. Int. J. Clim.; 2007; 27, pp. 1547-1578. [DOI: https://dx.doi.org/10.1002/joc.1556]

59. Tangdamrongsub, N.; Steele-Dunne, S.C.; Gunter, B.C.; Ditmar, P.G.; Weerts, A.H. Data Assimilation of GRACE Terrestrial Water Storage Estimates into a Regional Hydrological Model of the Rhine River Basin. Hydrol. Earth Syst. Sci.; 2015; 19, pp. 2079-2100. [DOI: https://dx.doi.org/10.5194/hess-19-2079-2015]

60. Tian, S.; Renzullo, L.J.; van Dijk, A.I.J.M.; Tregoning, P.; Walker, J.P. Global Joint Assimilation of GRACE and SMOS for Improved Estimation of Root-Zone Soil Moisture and Vegetation Response. Hydrol. Earth Syst. Sci.; 2019; 23, pp. 1067-1081. [DOI: https://dx.doi.org/10.5194/hess-23-1067-2019]

61. Vishwakarma, B.D.; Zhang, J.; Sneeuw, N. Downscaling GRACE Total Water Storage Change Using Partial Least Squares Regression. Sci. Data; 2021; 8, 95. [DOI: https://dx.doi.org/10.1038/s41597-021-00862-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33772016]

62. Arshad, A.; Mirchi, A.; Vilcaez, J.; Umar Akbar, M.; Madani, K. Reconstructing High-Resolution Groundwater Level Data Using a Hybrid Random Forest Model to Quantify Distributed Groundwater Changes in the Indus Basin. J. Hydrol.; 2024; 628, 130535. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2023.130535]

63. Chen, L.; He, Q.; Liu, K.; Li, J.; Jing, C. Downscaling of GRACE-Derived Groundwater Storage Based on the Random Forest Model. Remote Sens.; 2019; 11, 2979. [DOI: https://dx.doi.org/10.3390/rs11242979]

64. Zuo, J.; Xu, J.; Chen, Y.; Li, W. Downscaling Simulation of Groundwater Storage in the Tarim River Basin in Northwest China Based on GRACE Data. Phys. Chem. Earth Parts ABC; 2021; 123, 103042. [DOI: https://dx.doi.org/10.1016/j.pce.2021.103042]

65. 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, pp. 5973-5987. [DOI: https://dx.doi.org/10.1029/2017JD027468]

66. Sahour, H.; Sultan, M.; Vazifedan, M.; Abdelmohsen, K.; Karki, S.; Yellich, J.A.; Gebremichael, E.; Alshehri, F.; Elbayoumi, T.M. Statistical Applications to Downscale GRACE-Derived Terrestrial Water Storage Data and to Fill Temporal Gaps. Remote Sens.; 2020; 12, 533. [DOI: https://dx.doi.org/10.3390/rs12030533]

67. Yin, W.; Zhang, G.; Liu, F.; Zhang, D.; Zhang, X.; Chen, S. Improving the Spatial Resolution of GRACE-Based Groundwater Storage Estimates Using a Machine Learning Algorithm and Hydrological Model. Hydrogeol. J.; 2022; 30, pp. 947-963. [DOI: https://dx.doi.org/10.1007/s10040-021-02447-4]

68. Li, F.; Zhang, G.; Xu, Y.J. Assessing Climate Change Impacts on Water Resources in the Songhua River Basin. Water; 2016; 8, 420. [DOI: https://dx.doi.org/10.3390/w8100420]

69. Song, X.; Song, S.; Sun, W.; Mu, X.; Wang, S.; Li, J.; Li, Y. Recent Changes in Extreme Precipitation and Drought over the Songhua River Basin, China, during 1960–2013. Atmos. Res.; 2015; 157, pp. 137-152. [DOI: https://dx.doi.org/10.1016/j.atmosres.2015.01.022]

70. Zhang, B.; Liu, C.; Zhang, Z.; Xiong, S.; Zhang, W.; Li, Z.; An, B.; Wang, S. Downscaling and Calibration Analysis of Precipitation Data in the Songhua River Basin Using the GWRK Model and Rain Gauges. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2024; 17, pp. 12842-12853. [DOI: https://dx.doi.org/10.1109/JSTARS.2024.3424349]

71. Li, F.; Zhang, G.; Xu, Y.J. Spatiotemporal Variability of Climate and Streamflow in the Songhua River Basin, Northeast China. J. Hydrol.; 2014; 514, pp. 53-64. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2014.04.010]

72. Khorrami, B.; Ali, S.; Abadi, L.H.; Jehanzaib, M. Spatio-Temporal Variations in Characteristics of Terrestrial Water Storage and Associated Drought over Different Geographic Regions of Türkiye. Earth Sci. Inform.; 2023; 16, pp. 717-731. [DOI: https://dx.doi.org/10.1007/s12145-022-00907-3]

73. Khorrami, B.; Gunduz, O. Evaluation of the Temporal Variations of Groundwater Storage and Its Interactions with Climatic Variables Using GRACE Data and Hydrological Models: A Study from Turkey. Hydrol. Process.; 2021; 35, e14076. [DOI: https://dx.doi.org/10.1002/hyp.14076]

74. Khorrami, B.; Ali, S.; Sahin, O.G.; Gunduz, O. MODEL-COUPLED GRACE-BASED Analysis of Hydrological Dynamics of Drying Lake Urmia and Its Basin. Hydrol. Process.; 2023; 37, e14893. [DOI: https://dx.doi.org/10.1002/hyp.14893]

75. Rodell, M.; Houser, P.R.; Jambor, U.; Gottschalck, J.; Mitchell, K.; Meng, C.-J.; Arsenault, K.; Cosgrove, B.; Radakovich, J.; Bosilovich, M. et al. The Global Land Data Assimilation System. Bull. Am. Meteorol. Soc.; 2004; 85, pp. 381-394. [DOI: https://dx.doi.org/10.1175/BAMS-85-3-381]

76. Pekel, J.-F.; Cottam, A.; Gorelick, N.; Belward, A.S. High-Resolution Mapping of Global Surface Water and Its Long-Term Changes. Nature; 2016; 540, pp. 418-422. [DOI: https://dx.doi.org/10.1038/nature20584]

77. Chen, H. Hydrological Processes and Wetlands Distribution under Climate Change in the Amur River Basin: Impacts and Predictions. Ph.D. Thesis; Aerospace Information Research Institute, Chinese Academy of Sciences: Beijing, China, 2020.

78. Zhang, D.; Zhang, W. Study of Distributed Hydrological Model with the Dynamic Integration of Infiltration Excess and Saturated Excess Water Yielding Mechanism. Dr. Nanjing Nanjing Univ.; 2006; 529, 190.

79. Liu, C.; Xu, C.; Zhang, Z.; Xiong, S.; Zhang, W.; Zhang, B.; Chen, H.; Xu, Y.; Wang, S. Modeling Hydrological Consequences of 21st-Century Climate and Land Use/Land Cover Changes in a Mid-High Latitude Watershed. Geosci. Front.; 2024; 15, 101819. [DOI: https://dx.doi.org/10.1016/j.gsf.2024.101819]

80. Liu, Y.; Zhang, W.; Zhang, Z. A Conceptual Data Model Coupling with Physically-Based Distributed Hydrological Models Based on Catchment Discretization Schemas. J. Hydrol.; 2015; 530, pp. 206-215. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2015.09.049]

81. Wang, S.; Xu, C.; Zhang, W.; Chen, H.; Zhang, B. Human-Induced Water Loss from Closed Inland Lakes: Hydrological Simulations in China’s Daihai Lake. J. Hydrol.; 2022; 607, 127552. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2022.127552]

82. Xu, C.; Zhang, W.; Wang, S.; Chen, H.; Azzam, A.; Zhang, B.; Xu, Y.; Nie, N. Spatiotemporal Green Water Dynamics and Their Responses to Variations of Climatic and Underlying Surface Factors: A Case Study in the Sanjiang Plain, China. J. Hydrol. Reg. Stud.; 2023; 45, 101303. [DOI: https://dx.doi.org/10.1016/j.ejrh.2022.101303]

83. Zhang, D.; Zhang, W. Distributed Hydrological Modeling Study with the Dynamic Water Yielding Mechanism and RS/GIS Techniques. Remote Sensing for Agriculture, Ecosystems, and Hydrology VIII; SPIE: Stockholm, Sweden, 2006; Volume 6359, pp. 340-351.

84. Khorrami, B.; Pirasteh, S.; Ali, S.; Sahin, O.G.; Vaheddoost, B. Statistical Downscaling of GRACE TWSA Estimates to a 1-Km Spatial Resolution for a Local-Scale Surveillance of Flooding Potential. J. Hydrol.; 2023; 624, 129929. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2023.129929]

85. Senay, G.B.; Budde, M.; Verdin, J.P.; Melesse, A.M. A Coupled Remote Sensing and Simplified Surface Energy Balance Approach to Estimate Actual Evapotranspiration from Irrigated Fields. Sensors; 2007; 7, pp. 979-1000. [DOI: https://dx.doi.org/10.3390/s7060979]

86. Li, Q.; Shi, G.; Shangguan, W.; Nourani, V.; Li, J.; Li, L.; Huang, F.; Zhang, Y.; Wang, C.; Wang, D. et al. A 1 Km Daily Soil Moisture Dataset over China Using in Situ Measurement and Machine Learning. Earth Syst. Sci. Data; 2022; 14, pp. 5267-5286. [DOI: https://dx.doi.org/10.5194/essd-14-5267-2022]

87. Shangguan, W.; Li, Q.; Shi, G. A 1-Km Daily Soil Moisture Dataset over China Based on in-Situ Measurement (2000–2020); National Tibetan Plateau Data Center [Data Set]: Guangzhoug, China, 2022.

88. Ding, Y.; Peng, S. Spatiotemporal Trends and Attribution of Drought across China from 1901–2100. Sustainability; 2020; 12, 477. [DOI: https://dx.doi.org/10.3390/su12020477]

89. Peng, S. Spatiotemporal Change and Trend Analysis of Potential Evapotranspiration over the Loess Plateau of China during 2011–2100. Agric. For. Meteorol.; 2017; 233, pp. 183-194. [DOI: https://dx.doi.org/10.1016/j.agrformet.2016.11.129]

90. 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, pp. 1931-1946. [DOI: https://dx.doi.org/10.5194/essd-11-1931-2019]

91. Peng, S. 1-Km Monthly Mean Temperature Dataset for China (1901–2021); A Big Earth Data Platform for Three Poles: Lanzhou, China, 2019; [DOI: https://dx.doi.org/10.11888/Meteoro.tpdc.270961]

92. Arshad, A.; Mirchi, A.; Taghvaeian, S.; AghaKouchak, A. Downscaled-GRACE Data Reveal Anthropogenic and Climate-Induced Water Storage Decline Across the Indus Basin. Water Resour. Res.; 2024; 60, e2023WR035882. [DOI: https://dx.doi.org/10.1029/2023WR035882]

93. Abdi, H. Partial Least Squares Regression and Projection on Latent Structure Regression (PLS Regression). WIREs Comput. Stat.; 2010; 2, pp. 97-106. [DOI: https://dx.doi.org/10.1002/wics.51]

94. Vinzi, V.E. Handbook of Partial Least Squares; Springer: Berlin/Heidelberg, Germany, 2010.

95. Woldesenbet, T.A.; Elagib, N.A.; Ribbe, L.; Heinrich, J. Hydrological Responses to Land Use/Cover Changes in the Source Region of the Upper Blue Nile Basin, Ethiopia. Sci. Total Environ.; 2017; 575, pp. 724-741. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2016.09.124]

96. Wang, F.; Shen, Y.; Chen, T.; Chen, Q.; Li, W. Improved Multichannel Singular Spectrum Analysis for Post-Processing GRACE Monthly Gravity Field Models. Geophys. J. Int.; 2020; 223, pp. 825-839. [DOI: https://dx.doi.org/10.1093/gji/ggaa339]

97. Ghil, M.; Allen, M.R.; Dettinger, M.D.; Ide, K.; Kondrashov, D.; Mann, M.E.; Robertson, A.W.; Saunders, A.; Tian, Y.; Varadi, F. et al. ADVANCED SPECTRAL METHODS FOR CLIMATIC TIME SERIES. Rev. Geophys.; 2002; 40. [DOI: https://dx.doi.org/10.1029/2000RG000092]

98. Walwer, D.; Calais, E.; Ghil, M. Data-adaptive Detection of Transient Deformation in Geodetic Networks. J. Geophys. Res. Solid Earth; 2016; 121, pp. 2129-2152. [DOI: https://dx.doi.org/10.1002/2015JB012424]

99. Ali, S.; Liu, D.; Fu, Q.; Cheema, M.J.M.; Pham, Q.B.; Rahaman, M.M.; Dang, T.D.; Anh, D.T. Improving the Resolution of GRACE Data for Spatio-Temporal Groundwater Storage Assessment. Remote Sens.; 2021; 13, 3513. [DOI: https://dx.doi.org/10.3390/rs13173513]

100. Breiman, L. Random Forests. Mach. Learn.; 2001; 45, pp. 5-32. [DOI: https://dx.doi.org/10.1023/A:1010933404324]

101. Rahaman, M.M.; Thakur, B.; Kalra, A.; Li, R.; Maheshwari, P. Estimating High-Resolution Groundwater Storage from GRACE: A Random Forest Approach. Environments; 2019; 6, 63. [DOI: https://dx.doi.org/10.3390/environments6060063]

102. Brunsdon, C.; Fotheringham, A.S.; Charlton, M.E. Geographically Weighted Regression: A Method for Exploring Spatial Nonstationarity. Geogr. Anal.; 1996; 28, pp. 281-298. [DOI: https://dx.doi.org/10.1111/j.1538-4632.1996.tb00936.x]

103. Foody, G.M. Geographical Weighting as a Further Refinement to Regression Modelling: An Example Focused on the NDVI–Rainfall Relationship. Remote Sens. Environ.; 2003; 88, pp. 283-293. [DOI: https://dx.doi.org/10.1016/j.rse.2003.08.004]

104. Arshad, A.; Mirchi, A.; Samimi, M.; Ahmad, B. Combining Downscaled-GRACE Data with SWAT to Improve the Estimation of Groundwater Storage and Depletion Variations in the Irrigated Indus Basin (IIB). Sci. Total Environ.; 2022; 838, 156044. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2022.156044]

105. Chen, L.; Zhang, Z.; Liu, C.; Xiong, S.; Zhang, W.; Gao, H.; Yi, Y. Incorporating Glacier Processes into Hydrological Simulations in the Headwaters of the Yangtze and Yellow Rivers. Sci. Total Environ.; 2024; 951, 175474. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2024.175474]

106. Xu, C.; Zhang, Z.; Fu, Z.; Xiong, S.; Chen, H.; Zhang, W.; Wang, S.; Zhang, D.; Lu, H.; Jiang, X. Impacts of Climatic Fluctuations and Vegetation Greening on Regional Hydrological Processes: A Case Study in the Xiaoxinganling Mountains–Sanjiang Plain Region, Northeastern China. Remote Sens.; 2024; 16, 2709. [DOI: https://dx.doi.org/10.3390/rs16152709]

107. Chen, H.; Zhang, W.; Gao, H.; Nie, N. Climate Change and Anthropogenic Impacts on Wetland and Agriculture in the Songnen and Sanjiang Plain, Northeast China. Remote Sens.; 2018; 10, 356. [DOI: https://dx.doi.org/10.3390/rs10030356]

108. Yan, X.; Zhang, B.; Yao, Y.; Yang, Y.; Li, J.; Ran, Q. GRACE and Land Surface Models Reveal Severe Drought in Eastern China in 2019. J. Hydrol.; 2021; 601, 126640. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.126640]

109. Satish Kumar, K.; AnandRaj, P.; Sreelatha, K.; Bisht, D.S.; Sridhar, V. Monthly and Seasonal Drought Characterization Using GRACE-Based Groundwater Drought Index and Its Link to Teleconnections across South Indian River Basins. Climate; 2021; 9, 56. [DOI: https://dx.doi.org/10.3390/cli9040056]

110. Hong, Z.; Han, Z.; Li, X.; Long, D.; Tang, G.; Wang, J. Generation of an Improved Precipitation Data Set from Multisource Information over the Tibetan Plateau. J. Hydrometeorol.; 2021; 22, pp. 1275-1295. [DOI: https://dx.doi.org/10.1175/JHM-D-20-0252.1]

111. Long, D.; Pan, Y.; Zhou, J.; Chen, Y.; Hou, X.; Hong, Y.; Scanlon, B.R.; Longuevergne, L. Global Analysis of Spatiotemporal Variability in Merged Total Water Storage Changes Using Multiple GRACE Products and Global Hydrological Models. Remote Sens. Environ.; 2017; 192, pp. 198-216. [DOI: https://dx.doi.org/10.1016/j.rse.2017.02.011]

112. Cao, Y.; Nan, Z.; Cheng, G. GRACE Gravity Satellite Observations of Terrestrial Water Storage Changes for Drought Characterization in the Arid Land of Northwestern China. Remote Sens.; 2015; 7, pp. 1021-1047. [DOI: https://dx.doi.org/10.3390/rs70101021]