Improving Spatial Resolution of GRACE-Derived Water Storage Changes Based on Geographically Weighted Regression Downscaled Model

Obtaining high-resolution products that can accurately estimate the spatiotemporal changes in regional water resources is essential for the rational utilization of water resources. The temporal gravity field model derived from the Gravity Recovery and Climate Experiment (GRACE) satellite can effectively monitor the regional terrestrial water storage (TWS) changes. Nevertheless, the main of GRACE products is their coarse spatial resolution, which limits their applicability to small-scale regions. Herein, we propose a geographically weighted regression downscaled model (GWRDM) to improve the spatial resolution of TWS from 0.5° to 0.1° and validate it against the downscaled results provided by the spatial global regression downscaled model (SGRDM) and temporal grid regression downscaled model (TGRDM). The results show that the downscaled products based on GWRDM outperform those based on SGRDM and TGRDM. The GWRDM not only effectively improves the spatial resolution of GRACE products but also maintains the intensity of the original signal. Similarly, the GWRDM is further used to downscale groundwater storage (GWS) and validated against the in situ observations. The downscaled GWS based on GWRDM agrees well with in situ observations. On the monthly scale, the average correlation coefficient (CC) for both is 0.53 and above 0.70 for some in situ observations. On the annual scale, 81.25% and 37.50% of the in situ observations have CC values larger than 0.60 and 0.90, respectively, which further verifies the reliability of GWRDM. Consequently, the GWRDM provides a practical algorithm for obtaining high-resolution water storage estimates, which is expected to provide a reference for small-scale regional water resources management.


I. INTRODUCTION
I N RECENT years, drastic climate change has caused abnormalities in the global water cycle, which is prone to natural disasters, such as floods and droughts [1], [2]. Terrestrial water is an essential part of the worldwide hydrosphere cycle [3], which mainly includes soil water, groundwater, snow water equivalent, and canopy water. It is of great significance to understand the law of its changes for the rational allocation of water resources. Traditionally, monitoring of water resources has relied on hydrological stations and groundwater observation wells to obtain water level changes regularly. It is difficult to profoundly understand the water resource changes on a large scale [4]. There is a severe lack of traditional hydrological observation data, especially in mountainous areas with widespread deserts and complex topography [5]. In addition, the global hydrological models and land surface models have also been commonly used in water resource analysis over the past several decades. However, the models based on physical processes are sophisticated, and their reliability is usually affected by uncertain factors, such as climate-forcing data, model parameters, and model structures [6].
With the development of satellite observation technology, remote sensing methods have brought effective opportunities for water resource estimation. Noteworthy, the Gravity Recovery and Climate Experiment (GRACE) satellite, jointly developed by the National Aeronautics and Space Administration (NASA) and the Deutsches Zentrum für Luft und Raumfahrt, was successfully launched in March 2002, which provided an unprecedented perspective for detecting global and regional surface mass migration [7], [8], [9]. The temporal gravity field model derived from the GRACE observation is widely utilized to analyze the terrestrial water and groundwater storage (GWS) changes in medium-and large-scale regions [10], [11], [12], [13]. For example, Chen et al. [14] studied the water storage changes in the nine major basins of the Chinese Mainland in This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ detail through the Tongji-Grace2018 product combined with the GLDAS hydrological model, rainfall, and temperature data. The comparison with the in situ measurements reveals that the inversion results of GRACE have considerable reliability. Zhang et al. [15] used the GRACE model and 617 wells to analyze the GWS change in the North China Plain before and after the middle route of the South-to-North Water Diversion Project, showing that GRACE has an excellent ability to monitor groundwater changes.
However, limited by the performance of GRACE satellites, the spatial resolution of their products is still relatively coarse. GRACE products can only achieve high accuracy (1-1.5 cm equivalent water height) by inverting water storage changes in medium-and large-scale areas (∼200 000 km 2 ) [8], [16]. Therefore, GRACE products still have significant limitations in analyzing terrestrial water storage (TWS) changes in small-scale areas. Tapley et al. [17] showed that increasing the order of the GRACE spherical harmonic coefficient (SHC) product can improve its spatial resolution, but the noise of higher order coefficients also intensifies. To improve the accuracy of the SHC product, Gaussian filtering and destriping filtering algorithms were proposed to attenuate the high-frequency noise and the "stripe" error in the north-south direction [18], [19]. However, the averaging effect of the noise reduction algorithm will reduce the spatial resolution of the SHC. In addition to SHC products, agencies have released Mass Concentration (Mascon) products. Although the newly released Mascon version 06 product has improved the setting resolution, it still cannot satisfy the study of TWS changes in small-scale areas. Therefore, it is necessary to use effective methods to improve the spatial resolution of GRACE products.
In general, there are two broad categories for improving the spatial resolution of remote sensing products: enhancing the observation accuracy of instruments and downscaling methods. The GRACE satellite ended its mission in October 2017, and its successor mission, GRACE Follow-On, was successfully launched in May 2018. Previous research studies showed that the GRACE Follow-On satellite observation does not significantly improve spatial resolution [20]. Thus, it is difficult to improve the spatial resolution of the model by enhancing the observation methods. Therefore, downscaling techniques have been widely used to obtain high-resolution products of hydrological data in recent years [21], [22], [23]. Many attempts have been made to downscale the GRACE products [24], [25], [26], [27].
Nevertheless, the current GRACE downscaling research work is usually based on the temporal grid regression downscaled model (TGRDM) and the spatial global regression downscaled model (SGRDM). TGRDM establishes the regression relationship between GRACE and various hydrological variables on each coarse-resolution grid. It is considered that the regression relationship in the rough grid is suitable for high-resolution hydrological variables to realize the downscaling of GRACE products. However, the regression relationship established based on the grid shows strong constraints, so the downscaled signal has an apparent discontinuity between the coarse grids. The downscaled result still retains the characteristics of coarseresolution products, which do not conform to the continuous change of hydrological signals [28]. SGRDM refers to establishing the regression relationship between monthly GRACE and various hydrological variables in the spatial dimension. The downscaled signal has good continuity in space, but the signal strength has a specific attenuation. The reason is that when downscaling models are established globally, the individual grids are equally weighted. The influence of local features on the downscaled results is not considered. Several local or moving window methods have been used to improve the global downscaling model [29], [30], alleviating signal attenuation to a certain extent. However, local or moving window regression is still a spatial global regression model.
It is worth noting that the nonstationarity/spatial heterogeneity between GRACE and auxiliary hydrological variables is not adequately considered when downscaling the coarse-resolution GRACE products using TGRDM and SGRDM. According to the first theorem of geography proposed by Tobler [31], it can be known that "Everything is related to everything else, but near things are more related than distant things." Hydrological variables are usually spatially heterogeneous, so the influence of the surrounding grid on the central grid needs to be considered when downscaling the GRACE products. For the nonstationarity between geographical variables, Brunsdon et al. [32] proposed a geographically weighted regression (GWR) algorithm, which considers the spatial nonstationarity of local changes between independent and dependent variables.
Scholars have introduced the GWR algorithm into the hydrological variable downscaling work. For example, Xu et al. [33] explored the spatial heterogeneity between precipitation, normalized difference vegetation index (NDVI), and DEM through GWR. The monthly-scale precipitation datasets in the eastern Qinghai-Tibet Plateau and the Tianshan Mountains were downscaled from 0.25°(∼25 km) to 1 km to evaluate the accuracy of the GWR. And ground observations were used to validate the downscaled precipitation datasets. The results showed that the high-resolution precipitation dataset obtained by GWR outperformed traditional downscaling algorithms and had higher accuracy than the original dataset. Li et al. [34] downscaled the Tropical Rainfall Measuring Mission products in the Lancang River basin from 2001 to 2015 by GWR kriging, calibrated the downscaled products using geographical ratio analysis, and finally obtained satisfactory results. In addition, Duan et al. [35] proposed a GWR-based LST downscaling algorithm with better accuracy than the global TsHARP algorithm. Therefore, the GWR algorithm is widely applied in downscaling research of hydrological variables, such as land surface temperature, soil moisture, rainfall, and wind speed [36], [37], [38].
To alleviate the limitations of SGRDM and TGRDM, this study proposes a geographically weighted regression downscaled model (GWRDM) to downscale GRACE products. Moreover, the consistency of regional water storage estimates before and after downscaling is analyzed by correlation coefficient (CC), the Nash-Sutcliffe efficiency coefficient (NSE), and root mean square error (RMSE). Finally, we verify the reliability of the downscaled GWS through ground-based measurements. This study provides a tremendous potential approach for generating GRACE products with more satisfactory spatial

A. Study Area
The Tarim Basin and the Qaidam Basin are located between 34°41 N-45°06 N and 73°29 E-99°19 E in northwest inland China. Its climate is mainly warm temperate and highland continental, with less rainfall all year round and drought as the main characteristic. The Tarim Basin, located in southern Xinjiang, is the largest inland basin in China, situated between the Tianshan, Kunlun, and Altun mountains, with a high topography to the west and a low topography to the east [39]. The Qaidam Basin is located in northwestern Qinghai Province, northeastern Tibetan Plateau, northwest to the Altun Mountains, southwest to the Kunlun Mountains, and northeast to the Qilian Mountains [40]. Because of the special geographical location, the vast territory of the northwest region, and the extensive distribution of high mountains and deserts, it is challenging to deploy traditional hydrological stations [41]. There is a scarcity of observation that comprehensively reflects the dynamic changes in water resources. Therefore, satellite remote sensing technology is expected to be an essential tool for analyzing the water resource changes in the region.

B. Datasets
This study takes the TWS derived from the CSR M06 product as the target variable, and the auxiliary hydrological variables include precipitation (PRE), day land surface temperature (LST_day), night land surface temperature (LST_night), NDVI, and enhanced vegetation index (EVI) from Global Precipitation Mission (GPM), ERA_Land, and Moderate Resolution Imaging Spectroradiometer (MODIS) products. Generally, PRE, LST_day, and LST_night are essential factors affecting water resource changes [42]. Additionally, NDVI and EVI can effectively reflect the conditions of water resource changes [39], [43]. Thus, the above-mentioned hydrological variables are applied to downscale the GRACE-derived TWS. The summary of variables information is shown in Table I. 1) GRACE Product: The three official institutions, the Center for Space Research (CSR), Jet Propulsion Laboratory, and GeoForschungsZentrum Potsdam, have independently derived temporal gravity field models using GRACE satellite observation, mainly including SHC and Mascon products. The SHC products require data postprocessing, such as C20 replacement, degree 1 corrections, and Glacial Isostatic Adjustment correction. Moreover, adequate filter techniques are used to reduce high-frequency noise and "stripe" errors in the north-south direction [19]. However, we can analyze the variation of regional TWS directly through the Mascon product, which has been subjected to a priori constraints or regularization [44], [45]. Furthermore, previous studies have revealed that Mascon products have higher signal-to-noise ratios than SHC products [44], [46]. Therefore, to reduce the uncertainty of the postprocessing process, this study will utilize the CSR-M06 product to analyze the regional TWS changes [45], [47]. The CSR-M06 product has a temporal resolution of one month and a spatial resolution of 0.25°. The CSR-M06 has deducted the average value from 2004 to 2009 to obtain the terrestrial water storage anomalies (TWSA). Noteworthy, subsequent analyses will uniformly describe it as GRACE-derived TWS/TWSA. Mission to detect large-and medium-sized precipitation but also improves the detection of light precipitation (< 0.5 mm/h) and snowfall with higher spatiotemporal resolution and accuracy [49]. The GPM IMERG dataset achieves Scholars utilized MODIS products to analyze land, ocean, and atmospheric changes from regional to global scales; LST data comes from MOD11C3 version 6 products, mainly LST_day and LST_night. The MOD13C2 version 6 product provides NDVI and EVI, which are derived primarily from atmospherically corrected reflectance in the red, near-infrared, and blue bands. NDVI maintains the continuity of NOAA's AVHRR-NDVI time-series records for historical and climate applications, whereas EVI minimizes canopy-soil variation and increases sensitivity to densely vegetated conditions. These two products more effectively describe the state of vegetation on a global scale. The spatiotemporal resolutions of LST_day, LST_night, NDVI, and EVI are one month and 0.05°. 5) Ground-based Measurements: The data from 16 groundwater wells in northwest China from January 2005 to December 2015 were collected from the "China Geological Environmental Monitoring Groundwater Level Yearbook" compiled by the Ministry of Land and Resources. The red circles in Fig. 1 show the spatial distribution of well locations. Some observations had singular values since the early measured data were observed and recorded by humans. The measured data need to be preprocessed so that the water level change can more accurately reflect the regional groundwater change. The groundwater level change is obtained by subtracting the groundwater depth from the elevation. The water level change of each well is deducted from the average value from January 2005 to December 2015 to get the groundwater level anomalies (GWLA).

A. Brief Description of the Three Downscaling Models
The modeling forms of TGRDM, SGRDM, and GWRDM downscaling methods are given in Fig. 2, where M 1 , M 2 , M 3 , …, M n denote the spatial distribution of monthly independent or dependent variables within the study area, and the size of the monthly grid is h rows and k columns. TGRDM establishes the regression relationship between TWSA and hydrological variables in the temporal dimension for a certain grid [e.g., the red grid shown in Fig. 2(a)] and then gradually builds h × k downscaling models in the study area. The signals before and after downscaling in this way have good consistency in spatial distribution. Still, the downscaled signals are divided between adjacent coarse grids, which is inconsistent with the characteristics of continuous changes in hydrological signals. While SGRDM models the downscaling of GRACE products in the spatial dimension. This method establishes the regression relationship between TWSA and hydrological variables for a particular month (e.g., M 1 ) and then progressively builds downscaled models for n months [see Fig. 2(b)]. The downscaled signal based on SGRDM has a better continuity, but its amplitude has a more significant attenuation than the original signal. Previous studies have shown that both TGRDM and SGRDM have limitations due to the lack of consideration of the spatial heterogeneity in hydrological signals. Similarly, GWRDM performs downscaling modeling in the spatial dimension but takes into account the influence of the surrounding grid signal on the central grid signal to obtain downscaled results with higher accuracy. As shown in Fig. 2(c), the weight of the central grid signal is 1, and the further away from the central grid, the smaller the weight (that is, the darker the color, the larger the weight; conversely, the lighter the color, the smaller the weight, until it decreases to 0).

B. Construction of GWRDM
GWRDM captures the spatially nonstationary relationship between dependent and predictor variables by combining geographic coordinate information with linear regression methods [32] where y i denotes the dependent variable; x ik denotes the kth independent variable; (μ i , υ i ) denotes the geographic coordinates of the ith point; β 0 (μ i , υ i ) denotes the intercept of the ith point; β k (μ i , υ i ) denotes the coefficient of x ik ; and ε i denotes the residual of the ith point.
The regression parameter values in GWRDM vary with spatial location, and the parameters can be estimated by weighted least squares [32], [50] where X represents the independent variable matrix, w(u i , v i ) represents the weight matrix of the ith point, and other parameters are consistent with (1). The bisquare is used as the spatial kernel function to calculate the weight matrix of the model [32], which is shown as follows: where w ij represents the weight of observation point j, which is used to estimate the coefficient of observation point i; d ij represents the distance between observation points i and j; and b represents the adaptive bandwidth. When w ij =1, the GWR model degenerates into a global spatial model. The selection of an appropriate bandwidth parameter is crucial for the accuracy of the GWR estimation. Inappropriate bandwidth parameters will impact the accuracy of GWRDM results. Hence, the optimal modeling bandwidth is selected for the GWR model by the corrected Akaike Information Criterion (AICc) [51].
where n represents the number of sample points, S represents the hat matrix, tr(S) represents the trace of the hat matrix, and σ represents the estimate of the variance of the random error.
The bandwidth with the smallest AICc value is selected as the optimal bandwidth when establishing a monthly GWR model.

C. GWS Changes Estimated From GRACE and GLDAS
According to the water storage equations [52], [53], [54], the GRACE-derived TWSA is the sum of groundwater, soil moisture, snow water equivalent, and canopy water components. Thus, the groundwater storage anomalies (GWSA) can be obtained by deducting soil moisture, snow water equivalent, and canopy water from the TWSA as follows: (6) where SMSA represents the soil moisture storage anomalies, SWESA represents the snow water equivalent storage anomalies, and CWSA represents the canopy water storage anomalies. The sum of SMSA, SWESA, and CWSA is surface water storage anomalies (SWSA) obtained from the GLDAS Noah v2.1 model.

D. Evaluation Metrics
Three evaluation metrics are used to assess the performance of the above-mentioned downscaling schemes quantitatively. These include the CC, NSE, and RMSE, which are given as follows [55], [56], [57], [58]: where Y is the observed value, X is the predicted value,Ȳ andX are the mean of Y and X, respectively, and n is the number of datasets. For model building, the closer the CC and NSE between observed and predicted values are to 1, the better the model's accuracy. The smaller the RMSE, the closer the predicted value is to the observed value, and the higher the model's accuracy.

E. Research Framework
This research framework mainly includes data preprocessing and establishment of GWRDM, acquisition of high-resolution TWSA/GWSA downscaled products, and verification of the downscaled results, as shown in Fig. 3.

1) Data Preprocessing and Establishment of GWRDM:
To unify the spatial resolution of the study data, we aggregate the PRE, LST_day, LST_night, NDVI, and EVI to 0.1°a nd 0.5°, and the TWSA to 0.5°. Then, surface water is deducted from the TWSA to obtain 0.5°of GWSA. The GWRDM of TWSA/GWSA is established by combining the above-mentioned data with the 0.5°latitude and longitude information, and the 0.5°GWR coefficients and regression residuals of each covariate are available.

2) Acquisition of High-Resolution TWSA/GWSA Downscaled Products:
The scale is assumed to be invariant, meaning that the GWR relationships established at lowresolution scales still apply to other spatial resolutions [59]. Therefore, we employ the cubic spatial interpolation algorithm to transform GWR coefficients and residuals from 0.5°to 0.1°. The GWR coefficient of 0.1°and hydrological variables of 0.1°are used to calculate the TWSA/GWSA of 0.1°by substituting them into the GWR equation. Finally, the TWSA/GWSA downscaled results are obtained through regression residual correction.

3) Comparison and validation:
Meanwhile, the downscaled results of TGRDM and SGRDM are compared to those of GWRDM. Furthermore, we validate the downscaled results with the in situ observations.

IV. RESULT
A. Parameter Analysis of GWRDM 1) Bandwidth and Weight Distribution: Fig. 4 illustrates the monthly bandwidth and spatial distribution of weights for GWRDM. The differences between the bandwidths from month to month indicate that the relationship between TWSA and hydrological variables changes over time [see Fig. 4(a)]. The monthly bandwidth fluctuated wildly from 2013 to 2016, which may be caused by the rapid spatial changes of TWSA and hydrological signals. The enormous bandwidth is 69, and the smallest bandwidth is 54, accounting for the most significant proportion, about 85.10%. Fig. 4(b) presents the spatial distribution of weights for a grid (79.75°E-38.75°N) with a bandwidth of 54. It reveals that the farther from the center grid, the lower the importance. The sum of the weights within a 5 × 5 coarse grid (∼250 km × 250 km) accounts for about 81.76% of the total weights (purple dashed box). Similarly, the sum of the weights in the 7 × 7 coarse grid (∼350 km × 350 km) accounts for about 97.95% of the total weights (green dotted box), which agrees with the original resolution of the GRACE product [see Fig. 4(b)]. 2) Regression Coefficients: Fig. 5 presents the spatial distribution of the regression coefficient trends for each hydrological variable when building the GWRDM [see (1)]. The spatial distribution of each variable's GWR regression coefficients is quite different, indicating that the variables have different effects on the downscaling of TWSA. Furthermore, the regression coefficients of the same variable also have significant differences in the spatial distribution, demonstrating that the TWSA signal has spatial heterogeneity in this region.    mm/year, indicating that it effectively preserves the intensity of the original GRACE. Nevertheless, the signal jumps between adjacent grids and still has the grid characteristics of the coarse-resolution signal, which is not consistent with the continuous variation of TWS [see Fig. 6(j)-(l)]. The reason is that TGRDM calculates high-resolution downscaled signals through the same regression relationships within the same coarse-resolution grid. Meanwhile, the regression relationships between adjacent coarse-resolution grids are different. In contrast, the GWRDM downscaled results, which range from −46.112 to 34.982 mm/year, with an average of around −2.348 mm/year, preserve the spatial distribution of the GRACE well and have good continuity between adjacent coarse-resolution grids. Therefore, GWRDM, which considers the spatially heterogeneous relationship between explanatory variables and TWSA, produces optimal local downscaled results and effectively improves the deficiencies of SGRDM and TGRDM.

1) Spatial Distribution of Downscaled
2) Statistically Analysis of Accuracy: Due to the lack of high-resolution products to validate the downscaled results, we analyze the CC, NSE, and RMSE of the downscaled results' resampled signals and the GRACE in the spatiotemporal dimension. Fig. 7 presents the spatial CC, NSE, and RMSE values of GWRDM, SGRDM, and TGRDM. Meanwhile, Table II presents the proportion of grids with different accuracies. The CC values of SGRDM and TGRDM in the center were relatively small, whereas those in the north and south were above 0.8. The CC values of SGRDM and TGRDM between 0.8 and 1 were 67.04% and 91.12%, respectively. The performance of SGRDM on NSE is poor; the proportion of grids less than 0.6 is 46.81%, which is worse than that of GWRDM and TGRDM. The RMSE values of the SGRDM in the northwest, central, and southeastern regions were relatively large, which was related to the drastic signal changes. It indicates that the SGRDM could not sufficiently capture the rapidly changing signal characteristics of the local regions. The RMSE spatial distribution of TGRDM is similar to that of SGRDM, but its accuracy is better than that of SGRDM. The RMSE values of SGRDM and TGRDM between 0 and 40 mm accounted for 88.89% and 94.32%, respectively. The CC values of GWRDM are above 0.8, and the proportion in the range of 0.9-1 reaches 99.75%. Similarly, the NSE of GWRDM also has an outstanding performance, and the ratio in the field of 0.9-1 reaches 91.23%. Except for the middle of the study area, the grids with RMSE values of 0 to 20 mm accounted for 95.79%. Compared with SGRDM and TGRDM, the ratio of CC values in GWRDM from 0.8 to 1 increased by 32.96% and 8.88%, respectively; the ratio of NSE values at 0.8-1 increased by 77.18% and 59.9%, respectively; and the ratio of RMSE values at 0-20 mm increased by 39.05% and 26.62%, respectively.
To further compare the accuracy of the three methods, we flatten the monthly resampling signals and GRACE (that is, convert the h × k dimensional data to (h * k)×1 dimensional data) and analyze the CC, NSE, and RMSE of both (see Fig. 8).
Meanwhile, Fig. 8 also counts the maximum, minimum, mean, and standard deviation of CC, NSE, and RMSE. Overall, the CC of GWRDM varies from 0.96 to 0.99, with an average of around 0.98 and a standard deviation of around 0.01, outperforming SGRDM and TGRDM. The CC of SGRDM fluctuates around 0.77 but is slightly less stable than GWRDM, with a standard deviation of 0.06. The CC of TGRDM is between 0.36 and 0.97, with a standard deviation of 0.09, but it fluctuated wildly from 2005 to 2011. The NSE of the three methods has a similar performance to CC. The RMSE of GWRDM is the best, ranging from 3.03 to 18.17 mm, with a standard deviation of 3.46 mm. The RMSE performance of TGRDM is second, ranging from 8.84 to 43.87 mm, with a standard deviation of 7.01 mm. The RMSE of SGRDM is the worst, ranging from 9.07 to 60.04 mm, with a standard deviation of 11.09 mm. The analysis in the spatiotemporal dimension shows that the accuracy of GWRDM is better than that of SGRDM and TGRDM.

3) Effect of Regression Residual Correction: This section will
analyze the accuracy improvement of TWSA trends before and after residual correction. Fig. 9 presents the pixel scatter density plot of the trend between the GRACE [see Fig. 6(c)] and the downscaled result resampling signal. Fig. 9(a)-(c) shows the result without residual correction, and Fig. 9(d)-(f) shows the result with residual correction. After the residual correction, the accuracy of the three downscaled methods has been improved to a certain extent. The SGRDM has the most apparent improvement   Before the residual correction, the TGRDM performs the best because the grid signal of the TGRDM has a strong constraining effect on the downscaled model, which is also a fundamental reason for the discontinuity of the spatial distribution downscaled signal. After residual correction, the accuracy of GWRDM is still better than that of SGRDM and TGRDM. Thus, we utilize the downscaled results with the residual correction procedure for the subsequent analysis and discussion.

A. Comparison of Subregional Downscaled TWSA
The three downscaling models can achieve satisfactory accuracy when analyzing the TWSA time series on a basin scale (see Section IV-B). Swenson et al. [60] indicated that the larger the area of the study region, the higher the accuracy of the GRACE-derived TWS. According to the analysis presented in Section IV-A, the sum of the weights within a 5 × 5 (7 × 7) coarse grid accounts for about 81.76% (97.95%) of the total weights. Therefore, the study area is divided into subareas according to 5 × 5 and 7 × 7 coarse grids (Figs. 15 and 16) to analyze the downscaling accuracy of the three methods in small-scale areas.
Among them, 66 and 32 subregions of the two segmentation methods are selected for analysis (S1-S66 and S1-S32). Figs. 11 and 12 summarize the CC, NSE, and RMSE of the TWSA time series before and after downscaling in each subregion. The accuracy of GWRDM in each subregion is relatively stable, the RMSE value is close to 0 mm, and CC and NSE values are close to 1 and 1, respectively. The reason is that GWRDM is a local modeling algorithm that considers the spatial heterogeneity of TWSA and hydrological variables.
The CC and NSE of TGRDM are above 0.95 and 0.8, respectively. Moreover, the RMSE is less than 20 mm except for individual subregions. Similarly, TGRDM is a local modeling algorithm in theory. However, it does not consider the influence of the hydrological variable signals of the surrounding grids on the TWSA signals of the central grid, and its accuracy is slightly worse than that of GWRDM. The accuracy of SGRDM in subregions is the worst, and the CC/NSE of different subregions is quite different. In the 5 × 5 grid division, the RMSE of SGRDM is more than 60 mm in the area where the signal changes drastically (S27 and S28 in Fig. 15). It indicates that the downscaling model of TWSA established by SGRDM in large-scale regions cannot effectively recover the intense signal changes in local areas. Comparing Figs. 11 and 12 shows that the accuracy of GWRDM and TGRDM does not improve significantly with the increase in the subregion area. In contrast, the accuracy of SGRDM has shown a specific improvement. Therefore, we apply GWRDM for the subsequent discussion.
B. Applying GWRDM to GWSA 1) Analysis of Downscaled GWSA: To further verify the reliability of the downscaling method, we downscale the GWSA using a research framework similar to TWSA. Fig. 13 shows the spatial distribution of the GWSA trend from January 2003 to December 2016, before and after downscaling. We find that the groundwater reserves in the Tianshan region suffered a more significant deficit, whereas those in the southeastern Kunlun Mountains and the Qaidam Basin increased. Comparing Figs. 6(c) and 13(a) illustrates that the spatial distribution of GWS and TWS is consistent, revealing that GWS changes mainly because TWS changes in the study area. Fig. 13(b) shows the downscaled GWSA result with GWRDM, which shows that the spatial distribution of GWSA changes before and after downscaling is relatively consistent. Additionally, the downscaled GWSA shows a more detailed spatial signal. Fig. 13(c) illustrates the difference between Fig. 13(a) and the resampled signals of Fig. 13(b), which is mainly located in the region where the signal intensely varies. In Fig. 13(d)-(f),the inner coincidence accuracy of the downscaled resampled signal and the GRACE-derived GWSA is spatially analyzed (refer to Section IV-C). After the residual correction, the CC and NSE values of each grid are above 0.95 and 0.85, respectively, and the RMSE value of 99.54% of the grids is below 10 mm. Consequently, using GWRDM to downscale regional GWSA can also achieve reasonable accuracy.  2) Validation of Downscaled GWSA: Furthermore, the downscaled results are validated and analyzed using measured groundwater data. Data from 16 groundwater wells (shown as the red dot in Fig. 1) in Northwest China are collected to verify the downscaled GWSA results from 2005 to 2015. Note that to avoid introducing possible errors associated with unknown specific yields, we only analyze the measured data and downscaled results in terms of trends and CCs. Fig. 14 shows the GWLA of some groundwater wells and the downscaled GWSA of the grid in which they are located, which are in good agreement. Groundwater reserves in the northern Tarim Basin show a decreasing trend due to global warming factors and the retreat of the snow line in the Tianshan Mountains, which reduces groundwater recharge to the nearby areas.
In addition, groundwater abstraction is also an important reason [39]. The water levels of #15 and #16 located in the Qaidam Basin show an increasing trend, consistent with the downscaled GWSA trend in the area shown in Fig. 13(b). Previous studies have shown that changes in water resources in the Qaidam Basin are mainly affected by annual rainfall, evapotranspiration, and teleconnection factors (i.e., the Pacific Interdecadal Oscillation) [40], [61]. In particular, increased rainfall is an essential factor in increasing water storage. Fig. 14(f) counts the CC between the downscaled GWSA and the measured groundwater level change for 16 wells, including monthly and annual scales. The CC of all groundwater wells is positively correlated on the monthly scale and ranges from 0.265 to 0.769, with a mean value of 0.53. Among them, seven wells have a CC above 0.625 (#6, #9, #10, #11, #12, #13, and #14). Compared with the monthly scale, the CC of the annual scale increased. Except for #1, #5, and #15, the CC of groundwater wells is more significant than 0.60, with six wells having a CC above 0.90. The comprehensive analysis mentioned above indicates that the GWRDM established in this study can effectively improve the spatial resolution of TWS and GWS.

VI. CONCLUSION
The successful realization of the GRACE satellite mission provides an essential approach to studying the water resource changes in large-scale areas. However, the spatial resolution of its products is coarse and still has a significant limitation in the analysis of water resource changes in small-scale areas. Therefore, this study proposes the GWRDM to improve the spatial resolution of GRACE products to satisfy the research on small-scale regional water resource changes. We analyze in detail the accuracy of downscaled results estimated from GWRDM, SGRDM, and TGRDM. The primary conclusions are as follows.
1) The GWRDM of GRACE products was constructed in the spatial dimension by introducing a GWR algorithm to obtain the downscaled results of TWSA/GWSA with higher spatial resolution. The comparison of three downscaling algorithms shows that GWRDM effectively alleviates the discontinuity of downscaled signals between coarse grids compared with TGRDM. Meanwhile, GWRDM better preserves the intensity of the original signals than SGRDM. 2) In terms of spatial signal analysis, the CC, NSE, and RMSE of the downscaled results of GWRDM are better than those of SGRDM and TGRDM. The ratio of CC in GWRDM from 0.8 to 1 increased by 32.96% and 8.88%, respectively; the ratio of NSE at 0.8-1 increased by 77.18% and 59.9%, respectively; the ratio of RMSE at 0-20 mm increased by 39.05% and 26.62%, respectively. Furthermore, the accuracy of the downscaled results after residual correction is improved compared to before. Similarly, GWRDM has excellent accuracy in time series analysis and small-scale regional signal analysis.
3) The GWRDM is applied to downscale GWSA, and the results indicate that the signals before and after downscaling have good consistency in the spatial distribution. After residual correction, the spatial comparison of the downscaled signal with the original signal shows that the CC and NSE of each grid are above 0.95 and 0.85, respectively, and 99.54% of the grids have an RMSE of less than 10 mm. The in situ observations verified the downscaled GWSA. The results showed that the various trends are in good agreement, and the CC values of each measurement dataset show a positive correlation. On the annual scale, the CC of 81.25% of measurement data is more significant than 0.60. The CC of 37.50% of measurement data is more significant than 0.90, which further verifies that the downscaled results based on GWRDM are reliable. The effect of spatial dimensional weighting has been considered in the GWRDM, but there is no discussion of the effect on the temporal dimension. In future studies, the impacts of temporal dimensional weighting on the downscaling accuracy of GRACE products can be further analyzed.