SYSU TWSA v1.0: global high-resolution terrestrial water storage anomalies via satellite gravimetry

Publicly available global high-resolution terrestrial water storage anomaly (TWSA) datasets derived from satellite gravimetry remain scarce. Many existing global downscaling products rely heavily on hydrological models. Consequently, their performance can degrade in regions where key mass variations observed by the Gravity Recovery and Climate Experiment (GRACE) and its successor mission GRACE Follow-On (GFO) are poorly represented in the models, notably those associated with mountain glaciers and large lakes. Here we provide SYSU TWSA, a global monthly 0.5° TWSA dataset spanning April 2002 to December 2022, generated using a joint-inversion spatial downscaling framework that integrates large-scale constraints from GRACE/GFO, high-resolution spatial patterns from the WaterGAP Global Hydrological Model (WGHM), and additional mascon groups that explicitly represent mountain glaciers and selected large or rapidly changing lakes. The dataset helps alleviate the current shortage of global high-resolution products and explicitly strengthens the representation of glacier- and lake-related signals. We assess SYSU TWSA through four complementary evaluations: (1) basin-wise consistency with raw GRACE/GFO estimates, (2) a basin water-balance consistency check, (3) an independent evaluation against in situ groundwater well observations, and (4) comparisons with representative downscaled products in both the spectral and spatial domains. SYSU TWSA shows generally good agreement with GRACE/GFO at the basin scale, with coefficients of determination (R²) exceeding 0.85 across basin-size classes. In small basins, consistency with terrestrial water fluxes derived from the basin water-balance equation improves substantially, with NSE increasing by 17.1 % relative to raw GRACE/GFO across 1200 basins. Agreement with groundwater wells also improves, with correlations increasing at 67.7 % of 28 248 wells. Comparisons with representative assimilation-based and deep-learning downscaled products further indicate that SYSU TWSA demonstrates competitive overall accuracy while strengthening the representation of glacier- and lake-related signals. The SYSU TWSA dataset is openly available at the National Tibetan Plateau Data Center (https://doi.org/10.11888/Terre.tpdc.303322, Xiong et al., 2026).

1 Introduction

Water is a critical resource that sustains Earth's ecosystems and supports human societies, and its spatiotemporal distribution and availability are largely governed by the global water cycle. As an essential component of the climate system, the global water cycle redistributes water resources across space and time, thereby shaping regional water sustainability, the frequency and intensity of hydrological extremes (e.g., floods and droughts), and influencing global sea level change (Bierkens, 2015; Rodell et al., 2018; Reager et al., 2016). In this context, terrestrial water storage integrates changes across multiple water compartments, including surface water, soil moisture, groundwater, snow, and ice, and serves as a key indicator of regional water balance as well as the hydrological response to climate change and human activities. Consequently, high-accuracy, continuous monitoring of water storage dynamics is essential for advancing our understanding of the water cycle, supporting water-security management, and addressing global change.

With the launch of the Gravity Recovery and Climate Experiment (GRACE) mission and its successor, GRACE Follow-On (GFO), satellite gravimetry has provided continuous observations of time-variable gravity fields, enabling the tracking of global surface mass redistribution (Wahr et al., 1998; Tapley et al., 2019; Feng et al., 2023; Rodell and Reager, 2023). Among these mass transport processes, terrestrial water storage anomaly (TWSA) constitutes one of the dominant contributors to time-variable gravity signals (Rodell and Reager, 2023). However, constraints imposed by satellite orbit design, payload instrument precision, and background model errors limit the effective spatial resolution of monthly GRACE/GFO solutions to ∼330 km (Vishwakarma et al., 2018; Tapley et al., 2019), which substantially restricts their utility for small basin-scale hydrological research and applications (Pail et al., 2015; Chen et al., 2022; Wiese et al., 2022). Improving the effective spatial resolution of GRACE/GFO-derived TWSA is therefore critical for extending its hydrological applicability.

One of the main challenges in GRACE/GFO spatial downscaling is the lack of independent, high-precision measurements of TWSA, which makes it difficult to both develop downscaling frameworks and rigorously validate their accuracy. Existing studies generally follow two main pathways. The first is assimilation-based downscaling, which incorporates the large-scale observational constraints from GRACE/GFO into a high-resolution model to improve spatial resolution (Eicker et al., 2014; Li et al., 2019a; Gerdener et al., 2023a; Gou and Soja, 2024). Such approaches have important advantages because they can update model states and related hydrological fluxes in a physically consistent way (Springer et al., 2026). However, they typically rely heavily on structure and storage representation of the underlying hydrological models. Many models still exhibit large uncertainties in time-varying signals from interannual to long-term scales (Scanlon et al., 2018) and often have deficiencies or missing modules for key processes such as groundwater dynamics, mountain glacier mass changes and large-lake storage variations. This dependence is particularly important because GRACE/GFO observations only provide vertically integrated TWSA, whereas data assimilation must distribute the GRACE/GFO-derived increments among individual water storage compartments. If these compartments or processes are not assigned to the model, the increments may not be adequately represented in the model state vector, the increments may not be assigned to the correct storage components, limiting downscaling performance in regions where such processes dominate the observed mass signal (Gou and Soja, 2024). The second pathway is statistical downscaling. It establishes statistical relationships between GRACE/GFO observations and finer-scale variables (e.g., precipitation, evapotranspiration, and temperature) to distribute large-scale signals onto finer spatial scales (Ning et al., 2014; Yin et al., 2018; Vishwakarma et al., 2021). Recent advances in machine learning have further improved the ability to capture nonlinear relationships (Seyoum and Milewski, 2017; Zuo et al., 2021; Zhang et al., 2021; Yin et al., 2022). Nevertheless, global statistical downscaling often faces challenges related to the generalizability of empirical statistical relationships under spatial heterogeneity, overfitting in machine-learning training, and the reliability of transferring relationships from GRACE/GFO-effective scales to finer local scales. When the training samples are drawn directly from gridded GRACE/GFO TWSA products, such as official mascon products provided at 0.5 or 0.25° (Loomis et al., 2019; Save et al., 2016; Wiese et al., 2016), although distributed on fine grids, their information content is substantially coarser than the nominal resolution (Watkins et al., 2015; Vishwakarma et al., 2021). As a result, relationships learned from such targets may be valid only at large scales, that is, within the effective GRACE/GFO resolution, while their applicability at finer local resolutions, such as those of the high-resolution predictor variables, remains insufficiently demonstrated.

Against this backdrop, publicly available global high-resolution TWSA products remain scarce, and existing approaches still face distinct limitations. To address these issues, we recently developed a joint-inversion spatial downscaling framework and validated it across China (Xiong et al., 2025a). Conceptually, this framework differs from both data assimilation and statistical downscaling approaches. Unlike data assimilation, it does not update hydrological model states or require the GRACE/GFO-derived increments to be distributed among predefined model storage compartments. This avoids the need to modify the model state vector, specify model-error and observation-error covariance structures, or develop additional process modules when key mass signals are missing from the model. Unlike statistical or machine-learning downscaling, it does not train empirical relationships between GRACE/GFO TWSA and high-resolution predictors, thereby reducing the dependence on training samples, predictor selection, and scale-transfer assumptions from the effective GRACE/GFO resolution to finer grids. Instead, the framework uses WGHM and external datasets to define spatial basis functions, while the temporal evolution of these basis functions is estimated directly from GRACE/GFO observations. This formulation allows additional basis functions, such as glacier and lake mascon groups, to be incorporated when important mass signals are not adequately represented in the hydrological model. In this way, missing or poorly represented signals can be introduced into the inversion through dedicated spatial basis functions, without requiring explicit process modules or additional state variables within the model. Building on this framework, we extend the approach to the global domain and develop a high-resolution SYSU TWSA product to alleviate the current scarcity of global high-resolution datasets and to better represent glacierized regions and large lakes. We evaluate its consistency with GRACE/GFO at the basin scale, particularly for basins comparable to the effective GRACE/GFO resolution, and conduct independent checks using a water-balance constraint and in situ groundwater well observations to assess its potential advantages. Finally, we compare the SYSU product with existing representative global downscaled products.

2 Data

2.1 GRACE/GFO data and processing

In this study, we used the monthly GRACE/GFO Level 2 gravity field products from April 2002 to December 2022, provided in the form of spherical harmonic coefficients (up to degree and order 90), obtained from the Center for Space Research at the University of Texas (CSR). In this record, missing months were retained without interpolation or gap filling. The post-processing of these coefficients included the following steps: (1) replacing the degree-1 term with monthly geocenter estimates computed by Landerer (2019), and substituting the C₂₀ and C₃₀ coefficients with satellite laser ranging (SLR) estimates (Sun et al., 2016; Loomis et al., 2020); (2) removing the mean gravity field over January 2004 to December 2009; (3) applying a GIA correction using the ICE-6G-D model (Peltier et al., 2018); (4) filtering stripe noise using the DDK3 filter (Kusche et al., 2009); (5) reducing the topography effects (Yang et al., 2022); and (6) synthesizing the spherical harmonics to derive surface mass changes on a 0.5° grid (Wahr et al., 1998). Notably, although the surface mass changes were represented on a 0.5° grid after synthesis, their effective spatial resolution remained approximately 330 km.

2.2 WaterGAP global hydrological model outputs

The WaterGAP Global Hydrology Model (WGHM version 2.2e; Müller Schmied et al., 2024) outputs used in this study covered 2002–2022 and were forced by the GSWP3-ERA5 climate dataset. WGHM provides gridded monthly water storage components at 0.5° resolution, simulated within a vertical and lateral water-balance framework, and represents many storage compartments, including lakes, wetlands, rivers, reservoirs, groundwater, soil moisture, and snow water equivalent. Notably, WGHM does not simulate mountain glacier mass changes and does not explicitly represent several large lakes, such as the North American Great Lakes, Lake Victoria in Africa, and rapidly expanding lakes in the Tibetan Plateau.

2.3 Mountain glacier catalogue

The Randolph Glacier Inventory (RGI) version 6.0 is a global glacier outline dataset designed to provide a static snapshot of glacier extents worldwide (Pfeffer et al., 2014). It provides polygonal boundaries for individual glaciers, compiled through the Global Land Ice Measurements from Space (GLIMS) initiative. In this study, we used RGI v6.0 to delineate global mountain glacier mascon groups.

2.4 Hydrometeorological data

In this study, monthly precipitation (P), evapotranspiration (E), and runoff (R) data were adopted from the ERA5-Land dataset. These variables were used in the water-balance equation to evaluate the potential advantages of the SYSU product relative to the raw GRACE/GFO observations at the global level-4 basins (Lehner and Grill, 2013). ERA5-Land is a high-resolution land reanalysis produced by the ECMWF ERA5 system, and the data are available from the Copernicus Climate Data Store (Copernicus Climate Change Service, 2019).

2.5 In situ ground wells data

The in situ groundwater well observations used in this study were obtained from the Global Groundwater Monitoring Network (GGMN), coordinated by the International Groundwater Resources Assessment Centre (IGRAC), and were used as independent data to evaluate the performance of the SYSU TWSA product. It should be noted that GGMN does not provide a comprehensive global coverage of groundwater observations; its spatial distribution and temporal continuity depend on the degree to which national and sub-national agencies openly share their data, resulting in substantial heterogeneity across regions. To ensure comparability, we pre-processed the groundwater level time series by aggregating the original measurements to monthly means and retained only wells with at least 60 months (5 years) of valid records during April 2002 to December 2022. These records were not required to be fully continuous, as globally distributed well records without missing months are relatively scarce, and imposing a strict continuity requirement would have substantially reduced the number of wells available for validation. Therefore, a limited number of missing months was allowed, and missing values were not interpolated. In addition, anomalous monthly values, such as obvious outliers or short-term disturbances potentially related to pumping tests, were removed. This screening yielded 28 248 wells for the validation analysis.

Figure 1 Flowchart of the joint inversion downscaling method. SWES denotes snow water equivalent storage, SMS is soil moisture storage, GWS represents groundwater storage, and RESS refers to reservoir storage.

3 Joint inversion downscaling framework

The key strategy of the joint inversion downscaling was to integrate the large-scale information from GRACE/GFO with high-resolution spatial information from the model and self-defined mascon groups (Xiong et al., 2025a). The overall workflow was summarized in Fig. 1. TWSA can be represented as the sum of multiple water storage compartments. In this framework, each compartment is further expressed as a linear combination of spatial basis functions (S) and their temporal evolution coefficients (A), as follows:

$$\begin{array}{}\text{(1)}& \mathrm{TWSA}=\sum _{i=\mathrm{1}}^n{\mathrm{WC}}_i=\sum _{i=\mathrm{1}}^n\sum _{j=\mathrm{1}}^{m_i}{A}_{i,j}\cdot S_{i,j}\end{array}$$

where WC_i denotes the ith water storage compartment, n is the total number of water storage compartments, including snow water equivalent storage (SWES), soil moisture storage (SMS), groundwater storage (GWS), surface water storage (e.g., reservoir, river, and wet), mountain glacier, and lake. S_i,j is the jth high-resolution spatial basis function of component i, A_i,j is its corresponding temporal evolution coefficient, and m_i is the number of retained spatial basis functions for component i.

The high-resolution spatial basis functions were comprised of (i) spatial patterns extracted from WGHM-simulated TWSA compartments (Sect. 3.1) and (ii) supplemental basis functions constructed for signals not explicitly simulated in WGHM by defining global mountain glaciers and selected lakes mascon groups (Sect. 3.2 and 3.3). For the temporal evolution, we first ensured consistency with the effective GRACE/GFO resolution by truncating all spatial basis functions to the same maximum spherical harmonic degree and order as the GRACE/GFO fields and applying the DDK3 filter. Based on these filtered basis functions, the observation equation was written as follows:

$$\begin{array}{}\text{(2)}& {\mathrm{TWSA}}_{\mathrm{GRACE}/\mathrm{GFO}}=\sum _{i=\mathrm{1}}^n\sum _{j=\mathrm{1}}^{m_i}{A}_{i,j}\cdot S_{i,j}^{\mathrm{f}}+\mathit{\epsilon }\end{array}$$

where ${\mathrm{S}}_{i,j}^{\mathrm{f}}$ is the jth filtered spatial basis function of the ith water storage compartment, ε denotes the residual term. This scale matching is essential for the subsequent least-squares fit, because the temporal coefficients can be robustly estimated only when the observations and basis functions share a consistent effective spatial resolution. We then fitted GRACE/GFO observations using these filtered spatial basis functions via a least-squares adjustment to estimate the temporal coefficients (A_i,j) associated with each basis function. These adjusted coefficients replaced the potentially unreliable temporal evolution simulated by the model, and provided temporal information for the self-defined mascon groups. Finally, the downscaled TWSA fields were reconstructed by combining the adjusted temporal evolution (A) with the unfiltered (high-resolution) spatial basis functions (S). The final downscaled TWSA is provided on the 0.5° grid of the spatial basis functions derived from WGHM and self-defined mascon groups.

3.1 Extracting WGHM-simulated spatial patterns

Independent component analysis (ICA) is a spatiotemporal decomposition method based on high-order statistical information (Stone, 2004). It decomposes mixed signals into a set of spatial patterns and their corresponding temporal evolution coefficients, thereby representing TWSA as a linear combination of spatial basis functions and temporal coefficients. WGHM provides gridded monthly mass change for multiple vertically integrated TWSA compartments, including snow water equivalent, reservoir storage, soil moisture, groundwater, wets, and rivers. We applied ICA separately to each compartment to extract spatial patterns, which were then used as spatial basis functions in the inversion. During the ICA decomposition, we mainly retained the spatial patterns whose corresponding temporal coefficients represent trend, seasonal, and interannual variations, as these modes capture the main temporal behaviors of monthly water storage changes. Retaining only these dominant modes also reduces the number of unknown parameters, avoids an over-parameterized solution, and helps improve the stability of the inversion.

3.2 Defining spatial basis functions for glacier mascon groups

The glacier mascon groups are user-defined regions of the Earth's surface designed to encompass concentrated glacierized areas exceeding 100 km² (Ciracì et al., 2020). Each mascon group is represented by multiple small blocks, whose centroids are defined on the 0.5° grid. These mascon groups do not necessarily coincide with the exact glacier outlines, but are designed to capture glacier mass changes detectable by GRACE/GFO (Jacob et al., 2012).

Building on the global glacier mascon groups of Jacob et al. (2012) and Ciracì et al. (2020), we refined the mascon groups using the glacier information provided by RGI v6.0. Global mountain glaciers were divided into 14 major regions: (1) Alaska; (2) Western Canada and the United States; (3) Canadian Archipelago; (4) Iceland; (5) Svalbard; (6) Scandinavia; (7) Russian Arctic; (8) North Asia; (9) Central Europe; (10) Caucasus and Middle East; (11) High Mountain Asia; (12) Low Latitudes; (13) Southern Andes; and (14) New Zealand. Each major region was further subdivided into several subregions, shown in different colors (Fig. A1). Within each mascon unit, the glacier-covered areas were coded as 1 and the non-glacier areas as 0. We then incorporated these self-defined glacier mascon groups as spatial basis functions for global mountain glaciers, thereby supplementing the missing glacier component in WGHM within the joint-inversion framework.

3.3 Building spatial basis functions for lake mascon groups

We defined three lake mascon groups representing (1) the Great Lakes (North American), (2) Lake Victoria (Africa), and (3) 50 representative lakes (Tibetan Plateau), as shown in Fig. A2. The lake dataset in the Tibetan Plateau was taken from Li et al. (2019b). Because these lakes are relatively small, we applied ICA to the gridded lake storage variations and retained the leading spatial patterns that capture the long-term trend and the seasonal cycle. These patterns were then used as spatial basis functions for the lake group in the Tibetan Plateau. In contrast, the Great Lakes and Lake Victoria were defined as individual mascon units, assigned a value of 1 inside the lake and 0 outside. We then incorporated these lake mascon groups as spatial basis functions to represent lake storage variations observed by GRACE/GFO yet not explicitly simulated in WGHM.

4 Validation

In this section, we evaluate the performance of the SYSU downscaled product from four perspectives: (1) basin-scale consistency with raw GRACE/GFO estimates after effective-resolution matching; (2) assessment of its potential advantages using a water-balance equation; (3) independent evaluation against in situ groundwater well observations; and (4) comparison with other downscaled products in both the spectral and spatial domains. Because GRACE/GFO observations are used as constraints in the inversion, comparisons with GRACE/GFO are treated as consistency checks rather than independent validation.

4.1 Compared with GRACE/GFO at the basin scale

Given the lack of independent high-precision TWSA measurements at the global scale, a direct validation of downscaled products at the grid scale (0.5°) is challenging. The only globally available TWSA observations are provided by GRACE/GFO. However, GRACE/GFO does not truly resolve variability at the 0.5° grid scale and therefore cannot serve as a reliable reference at this scale. By contrast, GRACE/GFO observations are generally reliable at the basin scale when basin sizes are comparable to the effective GRACE/GFO resolution (Boergens et al., 2022), which enables an indirect evaluation. It should be noted, however, that a direct basin-scale comparison between a downscaled product and GRACE/GFO is still inappropriate. GRACE/GFO measurements can suffer from signal leakage in some basins (Landerer and Swenson, 2012). Although a range of signal restoration approaches has been proposed (Vishwakarma et al., 2017; Xiong et al., 2025b; Long et al., 2015), their accuracy is difficult to assess across individual basins in the absence of ground truth. Under these circumstances, it is more appropriate to forward-process the downscaled product by expanding it into spherical harmonics to the same maximum degree and order as GRACE/GFO and applying the same smoothing filter before evaluation. This ensures that the two datasets have comparable spectral content, that is, an equivalent effective resolution.

Figure 2 Global distribution of 89 river basins ranked by ascending area. Humidity is characterized by the Aridity Index (AI). Basin boundaries are based on GRDC (2020) datasets (available at: https://mrb.grdc.bafg.de/, last access: 1 July 2026).

In this study, we selected 89 basins of different sizes from the GRDC dataset (GRDC, 2020). Based on basin area, we grouped them into large basins (40 basins; >500 000 km²), medium basins (36 basins; 100 000–500 000 km²), and small basins (13 basins; 40 000–100 000 km²), as shown in Fig. 2. In addition, global hydroclimatic conditions were stratified into arid (A: 0.03 < AI < 0.2), semi-arid (SA: 0.2 < AI < 0.5), sub-humid (SH: 0.5 < AI < 0.65), and humid (H: AI > 0.65) classes based on the aridity index (Trabucco and Zomer, 2019).

We expanded the SYSU product into spherical harmonics up to the same maximum degree and order as GRACE/GFO and applied the same spatial smoothing filter, yielding a filtered SYSU product with an effective spatial resolution consistent with GRACE/GFO. To assess performance across basin scales, we averaged (latitude-weighted) both GRACE/GFO and the filtered SYSU TWSA for large, medium, and small basins to derive basin-wise time series. We then concatenated the time series within each basin-size class into column vectors and generated scatter-density plots with GRACE/GFO on the X-axis and the filtered SYSU product on the y axis, from which the coefficient of determination (R²) was calculated. As shown in Fig. 3, after matching the effective spatial resolution, the filtered SYSU product agrees closely with GRACE/GFO in basin-mean TWSA, with R² exceeding 0.85 for all basin-size classes. Specifically, R² reached 0.95 for large basins, 0.87 for medium basins, and 0.86 for small basins. These results indicate that the SYSU product is highly reliable and support the effectiveness of the joint inversion downscaling framework in preserving the basin-scale accuracy of GRACE/GFO signals. Notably, performance is relatively lower for small basins. This is expected because GRACE/GFO observations are generally more reliable for large basins (Vishwakarma et al., 2018; Boergens et al., 2022), whereas small basins may perform relatively poorly due to their inherently low signal strength, which is more significantly affected by the measurement accuracy of GRACE/GFO (∼2 cm) and leakage errors.

Figure 3 Scatter density plots of the basin-wise average TWSA in different scale basins estimated from the filtered SYSU products (y axis) versus GRACE/GFO measurements (x axis).

As shown in Fig. 4, we compared the TWSA time series for different scale basins. For large basins, both SYSU and filtered SYSU products agree well with GRACE/GFO observations. The median root-mean-square error (RMSE) between the filtered SYSU product and GRACE/GFO is 1.27 cm, which falls within the GRACE/GFO measurement uncertainty (Vishwakarma et al., 2018; Tapley et al., 2019). The corresponding median Nash–Sutcliffe efficiency (NSE) is approximately 0.95, and the median correlation coefficient (CC) is 0.98. For medium basins, performance decreases slightly relative to large basins, but the majority of basins still show strong agreement in their time series, with an RMSE of 2.36 cm, an NSE of 0.83, and a CC of 0.95. For small basins, the consistency decreases further, with an RMSE of 3.75 cm, an NSE of 0.81, and a CC of 0.95.

Figure 4 Comparison of time series in different basin scales: (a) large basins, (b) medium basins, and (c) small basins.

Although GRACE/GFO observations are generally affected by signal attenuation and leakage, these effects are relatively small for large basins (Xiong et al., 2025b). Spherical harmonic truncation and spatial smoothing can cause signal spreading beyond basin boundaries, but most of the signal remains concentrated within the target region. As the joint-inversion downscaling framework preserves the large-scale accuracy of GRACE/GFO measurements, all three products show a high level of consistency in large basins. However, as the basin area decreases, the effects of attenuation and leakage typically become more pronounced due to the limited effective spatial resolution of GRACE/GFO. This tends to dampen basin-wise TWSA amplitudes and underestimate long-term trends (Scanlon et al., 2012; Long et al., 2015; Vishwakarma et al., 2016). Improving spatial resolution can alleviate these issues. Consequently, in basins where attenuation and leakage are severe, the basin-average TWSA derived from the SYSU product may exhibit larger amplitudes and trends than GRACE/GFO, as observed for the Cuyuni River basin (ID: 10), the Daugava Dvina basin (ID: 11), the Gota basin (ID: 3), and the Glomma basin (ID: 1). Consistently, when SYSU is degraded to the same effective resolution as GRACE/GFO (i.e., the filtered SYSU product), the differences in basin-average TWSA are substantially reduced. In addition, for basins in arid regions, such as the Hamun basin (ID: 40) and Lake Balkhash (ID: 45), the filtered SYSU product still shows discrepancies relative to GRACE/GFO. This may reflect the weak hydrological signals in these arid environments, for which both GRACE/GFO and WGHM have limited sensitivity to detect or simulate the variations accurately (Gou and Soja, 2024; Xiong et al., 2025a).

Figure 5 Comparison of TWSA signals across different basin scales. (a) Scatter plots comparing long-term trends, annual amplitudes, and semi-annual amplitudes estimated by the SYSU and filtered SYSU products (y axis) against GRACE/GFO observations (x axis). (b) Spatial distribution map illustrating basin-scale trend of the SYSU product and GRACE/FO observations. (c) Spatial distribution map showing basin-wise annual amplitude of the SYSU product and GRACE/FO observations.

To verify the performance of the SYSU product in preserving long-term trends and the amplitudes of the annual and semi-annual signals, we compared the consistency of signal estimates among the three TWSA products in large, medium, and small basins, i.e., the SYSU, the filtered SYSU, and GRACE/GFO, as shown in Fig. 5. For large basins, the long-term trend and the annual and semi-annual amplitudes estimated from SYSU closely match those from GRACE/GFO, with R² of 0.90, 0.99, and 0.96, respectively, and most points clustering near the red dashed line. Notably, the filtered SYSU product does not show a clear improvement over the SYSU in its agreement with GRACE/GFO, with R² of 0.94, 0.99, and 0.92, respectively. This indicates that GRACE/GFO signals are reliable for large-basin scales, where leakage effects are relatively small, and that our downscaling framework effectively inherits their accuracy. For medium basins, the agreement between SYSU and GRACE/GFO decreases, with R² of 0.69, 0.94, and 0.83 for the long-term trend, the annual and semi-annual amplitudes, respectively. The filtered SYSU product also shows no substantial improvement, with R² of 0.72, 0.97, and 0.80. In contrast, for small basins, SYSU shows relatively poor agreement with GRACE/GFO observations, whereas the consistency between the filtered SYSU product and GRACE/GFO increases markedly, with R² of 0.70, 0.90, and 0.95. This contrast indirectly suggests that GRACE/GFO performs poorly at small basins, where severe leakage errors can prevent the native observations from capturing the underlying signals. By increasing the spatial resolution, the SYSU product can substantially alleviate this problem. Therefore, the apparent mismatch between SYSU and GRACE/GFO, together with the closer agreement between filtered SYSU and GRACE/GFO, mainly reflects differences in effective spatial resolution.

Figure 5b and c compares the basin-wise secular trends and annual amplitudes derived from the SYSU product and GRACE/GFO. Overall, the two datasets agree closely, with RMSEs of 0.24 cm yr⁻¹ for the trend and 1.93 cm for the annual amplitude. The largest trend difference is found in the Daugava Dvina basin (ID: 11), reaching −0.67 cm yr⁻¹, whereas the largest discrepancy in annual amplitude occurs in the Cuyuni basin (ID: 10), at 15.23 cm. These differences are consistent with leakage effects (Fig. 4).

Figure 6 Agreement assessment of TWF estimates in 1200 level-4 basins. Comparison between SYSU product, GRACE/GFO observations, and water balance-derived TWF quantified via NSE. (a) Basin-wise NSE of SYSU product. (b) Basin-wise NSE of GRACE/GFO observations. (c) NSE differences between SYSU and GRACE/GFO across watersheds sorted by ascending drainage area; positive values indicate superior performance of the SYSU product. The red dashed line marks the GRACE/GFO effective resolution limit of 63 000 km² (Vishwakarma et al., 2018).

4.2 Evaluation based on water balance equation in sub-basins

To evaluate the potential advantages of the downscaled product, we used a water balance equation ($\mathrm{TWF}=P-E-R$) to assess whether the SYSU product improves performance relative to the raw GRACE/GFO observations. We selected 1200 level-4 basins globally (Lehner and Grill, 2013) as the evaluation domain and used NSE as the performance metric. Considering data gaps in the GRACE/GFO record and the degraded data quality during the single-accelerometer period, we restricted the analysis to January 2003–December 2010. For the SYSU and GRACE/GFO, terrestrial water flux (TWF) was estimated from monthly TWSA using a centered difference approach ($\mathrm{TWF}=\mathrm{TWSA}(t+\mathrm{1})-\mathrm{TWSA}(t-\mathrm{1})/\mathrm{2}\mathrm{\Delta }t$, (Long et al., 2014). For ERA5-Land, TWF was computed from P, E, and R, and an adjacent-month smoothing scheme was applied to reduce potential high-frequency noise (Landerer et al., 2010).

Figure 6 summarizes the agreement between TWF estimates derived from the SYSU product and GRACE/GFO, and those computed from an independent water balance equation. The SYSU product shows improved consistency with basin-wise water balance closure across the 1,200 basins, with 71.7 % of basins yielding positive NSE values. By comparison, 59.2 % of basins show positive NSE when TWF is derived from the raw GRACE/GFO. As shown in Fig. 6a and b, most negative NSE values occur in arid regions where hydrological signals are weak, such as the Sahara in Africa and the Arabian Peninsula in the Middle East. For basins larger than the GRACE/GFO effective-resolution limit (63 000 km²), SYSU provides little improvement over GRACE/GFO, and the two products yield nearly identical NSE values. In contrast, for basins smaller than this limit, the NSE increases by 17.1 %, indicating a potential advantage of SYSU for water balance consistency at finer spatial scales. The most evident improvements relative to GRACE/GFO are found in small basins in regions such as the western European coasts and the eastern Chinese coasts (Fig. 6a-b). These results suggest that, provided P, E, and R estimates are sufficiently accurate, improving the spatial resolution of TWSA products can yield TWF estimates that are closer to those derived by the basin water balance equation. It is worth noting that this comparison should be interpreted as an independent consistency check rather than a strict validation, because uncertainties in P, E, and R and unaccounted human interventions (e.g., inter-basin transfers) can affect basin water-balance closure.

Figure 7 Difference in correlation coefficients between the SYSU product and in situ groundwater well observations versus raw GRACE/GFO and well observations, calculated across 28,248 groundwater wells. Positive values indicate stronger agreement between the SYSU product and well data.

4.3 Comparison with in situ wells data

To evaluate the extent to which the SYSU TWSA product improves spatial detail relative to the coarse resolution GRACE/GFO observations, we used 28 248 in situ groundwater wells worldwide as an independent reference. Because well records measure groundwater level changes whereas satellite gravimetry provides water storage changes, a strict magnitude comparison would require converting water levels to storage using parameters such as specific yield. Such parameters are difficult to obtain reliably at the global scale. We therefore used the correlation coefficient as a consistency metric to quantify whether SYSU better matches well observations than the raw GRACE/GFO product. Specifically, for each well location, we interpolated the SYSU and GRACE/GFO TWSA fields to the site to obtain corresponding time series, and then computed their correlations with the observed groundwater level time series. The statistical significance of the correlations was assessed at the 95 % confidence level, and only significant correlations were retained for the subsequent comparison. By comparing correlations before and after downscaling, we tested whether the added spatial detail improves the consistency between TWSA estimates and local groundwater-level variations. As shown in Fig. 7, downscaling leads to a clear improvement: compared with raw GRACE/GFO, correlations between SYSU TWSA and groundwater level time series increase at 67.7 % of wells. Regional statistics indicate particularly pronounced improvements in the United States, India, and Australia, where correlations increase at 77.4 %, 76.2 %, and 77.7 % of wells, respectively, whereas the improvement is less evident in Europe. Overall, the results indicate that the downscaled SYSU TWSA is more consistent with observed groundwater level variations than the raw GRACE/GFO observations, although this correlation-based comparison should not be interpreted as a direct validation of groundwater storage changes.

4.4 Comparison with other products

We further compared the SYSU product with GRACE/GFO spherical harmonic and mascon solutions, as well as with WGHM and several representative downscaled products, to evaluate the performance of each dataset. The downscaled products considered here were derived from WGHM through data assimilation. One uses a conventional Ensemble Kalman Filter (EnKF) approach (GLWS 02, hereafter GLWS), and the other is based on a self-supervised deep learning approach (Gou and Soja, 2024). The evaluation was conducted from two perspectives, in the spectral and spatial domains.

Figure 8 RMS of spherical harmonic coefficients derived by different TWSA products (2003–2019). RMS is expressed in terms of geoid change here, in mm-units.

4.4.1 Spectral domain

In this section, we compare each downscaled product with GRACE/GFO and WGHM in the spherical harmonic domain. To ensure a consistent evaluation, we converted land TWSA from all downscaled products, the WGHM outputs, and GRACE/GFO (including both spherical harmonic and mascon solutions) into spherical harmonic coefficients. Notably, the expansion was performed over the global land domain and excluded the Greenland and Antarctic ice sheets to avoid contamination of the assessment by polar mass variations.

Figure 8 shows the RMS of each spherical harmonic coefficient's contribution to geoid variability (up to degree and order 90) for the different products. In general, in the spectral domain, lower degrees correspond to longer-wavelength, large-scale mass variability, whereas higher degrees correspond to shorter-wavelength variability. The results indicate that all downscaled products effectively remove much of the high-frequency content inherited from WGHM, yet their RMS values remain slightly higher than those of GRACE/GFO. In addition, compared with the EnKF-based GLWS product, both the deep learning downscaling product of Gou and Soja (2024) and the joint-inversion SYSU product better preserve the low-degree coefficients observed by GRACE/GFO. However, beyond roughly degree 30, SYSU retains substantially more variability than Gou and Soja's product.

Figure 9 shows the degree variance curves for October 2005, and the overall behavior is consistent with the patterns in Fig. 8. A clear offset is evident between WGHM and GRACE/GFO. This may reflect stronger signals in WGHM over Siberia and North America (Gerdener et al., 2023a), or may arise because the WGHM fields were not demeaned. Over degrees 0–25, all downscaled products remain in close agreement with GRACE/GFO. Beyond about degree 30, however, the curves begin to diverge. Because GRACE/GFO solutions are filtered and regularized, their high-degree power is suppressed. In contrast, the SYSU and GLWS products retain more high-frequency variance than the product of Gou and Soja (2024), indicating that they preserve more small-scale signals but also carry a higher noise level.

Figure 9 Degree variance for different products in October 2005 expressed by geoid heights (mm).

Figure 10 Time series of C₂₀ spherical harmonic coefficient derived by different products.

Given that individual spherical harmonic coefficients can be used to assess how a downscaling method corrects the hydrological signal from WGHM (Gerdener et al., 2023a), we selected C₂₀ as a representative example (Fig. 10). The GLWS product is closer to GRACE/GFO than WGHM, but it still tends to overestimate GRACE/GFO and exhibits larger temporal variability. This suggests that the EnKF-based assimilation inherits characteristics from both GRACE/GFO and WGHM and improves the hydrological model to some extent, yet remains constrained by its underlying error-model assumptions (Gerdener et al., 2023a). By contrast, the downscaled product of Gou and Soja (2024) and the SYSU product agree closely with GRACE/GFO, with C₂₀ values lying between the GRACE/GFO spherical harmonic and mascon solutions. These results indicate that deep learning-based and joint-inversion downscaling can offer advantages over conventional EnKF assimilation.

Figure 11 (a) Long-term trends in TWSA estimated from different products for the period January 2003 to December 2019. The inset maps highlight representative glacierized regions, including Alaska and the Tibetan Plateau (TP), where WGHM does not explicitly simulate mountain glacier mass changes. (b) Same as (a), but for annual amplitudes. (c) Same as (a), but for phases, indicating the month of peak amplitude (e.g., 1 represents January).

4.4.2 Spatial domain

In this section, we focus on spatial comparisons. Figure 11a shows the long-term trends estimated from each product from January 2003 to December 2019. Consistent with previous studies (Döll et al., 2014; Scanlon et al., 2018), the WGHM-simulated TWSA trends are substantially smaller in magnitude than those derived from GRACE/GFO. In addition, except for the polar regions, WGHM does not effectively capture the mountain glacier melt signals evident in GRACE/GFO observations. Despite these discrepancies, the two datasets show good agreement in some regions, particularly in areas with intensive groundwater depletion (e.g., the Northwest India Plain, the North China Plain, and California's Central Valley) and in regions characterized by pronounced wetting-drying variability (e.g., the Amazon and La Plata basins and the African Great Lakes region). Notably, relative to GRACE/GFO, WGHM exhibits richer spatial detail in the trend patterns.

All three downscaled products largely retain the spatial resolution of the WGHM. Their long-term trends exhibit much finer spatial detail than GRACE/GFO, revealing more localized positive and negative trend patterns. For example, the river-channel structure within the Amazon Basin is clearly resolved. However, the products differ in their ability to correct WGHM-related biases in regions where the model is structurally limited. In some regions, such as the northern High Plains of the United States, GRACE/GFO observations indicate a recovery in groundwater storage (Rateb et al., 2020; Scanlon et al., 2023), whereas the downscaled product of Gou and Soja (2024) does not correct the WGHM simulation and still suggests continued groundwater depletion. By contrast, the SYSU and GLWS products successfully adjust the negative WGHM trend and align with the GRACE/GFO observations, underscoring the contribution of GRACE/GFO constraints within both the joint-inversion framework and the EnKF approach.

More importantly, relative to the other two downscaled products, SYSU incorporates mascon groups for mountain glaciers and selected lakes during the inversion, enabling it to capture glacier-melt signals and storage increases in large lakes. This advantage is particularly evident in regions where WGHM lacks the corresponding storage components. In glacierized regions, such as Alaska and the Tibetan Plateau, WGHM does not explicitly simulate mountain glacier mass changes and therefore fails to represent the strong negative trends observed by GRACE/GFO. The GLWS and Gou and Soja (2024) products partly reflect these large-scale glacier-related signals, but their spatial patterns remain less localized or less complete in some glacierized regions. By contrast, SYSU explicitly incorporates glacier mascon basis functions and therefore provides a more localized representation of glacier mass loss while remaining broadly consistent with the large-scale GRACE/GFO signal.

A similar issue occurs for several large lakes that are detectable by GRACE/GFO but are not explicitly represented in WGHM. By incorporating dedicated lake mascon basis functions, SYSU better represents storage variations in selected large-lake regions, such as the North American Great Lakes and Lake Victoria. These results indicate that the supplemental glacier and lake basis functions help reduce signal distortion in regions where WGHM lacks the corresponding storage components. They also demonstrate that the joint-inversion framework can flexibly incorporate additional spatial basis functions for poorly represented mass-change processes, rather than relying solely on WGHM-derived spatial patterns. In contrast, in arid regions where hydrological signals are weak, such as northern Africa and the Middle East, differences among products are more difficult to interpret because both GRACE/GFO observations and hydrological-model simulations have larger relative uncertainties.

Figure 11b and c compares the annual amplitude and its phase (defined as the month of maximum amplitude) among the different products. Because seasonal water storage variation in most basins is dominated by the annual signals, the products show broadly consistent annual amplitudes. Moreover, the downscaled products not only achieve a spatial resolution comparable to that of WGHM, but also produce amplitudes more consistent with GRACE/GFO than WGHM. This indicates that these downscaling approaches benefit from incorporating GRACE/GFO constraints to improve hydrological-model representations of seasonal variation. For the phase (Fig. 11c), WGHM exhibits a 1–2 month offset relative to GRACE/GFO in multiple regions. GLWS reduces this discrepancy to some extent but still deviates from GRACE/GFO, whereas the SYSU and Gou and Soja (2024) products show phases that are generally closer to GRACE/GFO. Locally, however, phase advances or delays are still apparent, which may reflect genuine regional hydrological differences. A notable discrepancy appears in northern Africa, where the phase from SYSU and from Gou and Soja (2024) differs substantially. SYSU is close to the GRACE/GFO spherical harmonic solution, whereas Gou and Soja (2024) is close to the mascon solution. This difference likely reflects the different GRACE/GFO inputs used by the two products, with SYSU based on spherical harmonic coefficients and Gou and Soja (2024) based on the mascon product. The contrast may be further amplified in northern Africa, where extensive desert cover leads to weak hydrological signals and relatively large uncertainty in GRACE/GFO observations (Boergens et al., 2022). Moreover, given the limited effective spatial resolution of GRACE/GFO, it remains difficult to determine which downscaled product provides a more accurate phase estimate in this region.

Figure 12 Basin-wise RMSE comparison for three downscaled products and the WGHM output compared to GRACE/GFO observations.

Table 1 Evaluation of different products in 288 basins. Metrics include the median values of RMSE (in cm), NSE, CC, and R². For each row, the best-performing value is shown in bold, and the second-best is underlined.

The qualitative comparisons of long-term trends, annual amplitudes, and phases suggest that EnKF, deep learning, and our joint-inversion downscaling approaches all have considerable potential to improve WGHM-simulated TWSA estimates. Given that GRACE/GFO is highly accurate at its effective spatial resolution (Vishwakarma et al., 2018; Boergens et al., 2022), we quantify the improvement of each of the three downscaled products over WGHM in estimating TWSA. Specifically, we assess which product matches GRACE/GFO most closely when evaluated at the same effective resolution. To this end, we expanded the WGHM outputs and the three downscaled products into spherical harmonics to the same maximum degree and order as GRACE/GFO (degree/order 90) and applied the same DDK3 filtering to match the effective spatial resolution of GRACE/GFO. We then derived basin-average time series for 288 global level-3 basins (Lehner and Grill, 2013) over January 2003 to December 2019 and evaluated their agreement with GRACE/GFO using RMSE, NSE, CC, and R² (Table 1).

In terms of RMSE, all three downscaled products outperform WGHM (Table 1 and Fig. 12). Specifically, the median RMSE between WGHM and GRACE/GFO is 3.27 cm (mean: 4.41 cm). The GLWS product performs similarly to WGHM overall, with a median RMSE of 3.40 cm (mean: 4.04 cm). By contrast, the SYSU product and the Gou and Soja (2024) product achieve low RMSEs and both outperform GLWS. SYSU yields a median RMSE of 2.24 cm (mean: 2.44 cm), whereas Gou and Soja (2024) yields a median RMSE of 1.87 cm (mean: 2.37 cm). These values support the overall quality of the latter two downscaled products, given that GRACE/GFO uncertainties are on the order of 2–3 cm (Wahr et al., 2006; Tapley et al., 2019). Notably, in glacierized regions, SYSU performs markedly better than the Gou and Soja (2024) product (e.g., along the Gulf of Alaska coast). This advantage likely reflects that SYSU not only leverages information from WGHM but also explicitly incorporates mountain-glacier mascon groups. Finally, the largest RMSE values for SYSU (11.66 cm) and for Gou and Soja (2024) (17.06 cm) both occur in the glacierized Canadian Arctic Archipelago, which may be influenced by leakage from Greenland.

Figure 13 Same as Fig. 12, but for NSE.

Figure 14 Same as Fig. 12, but for CC.

In terms of NSE (Table 1 and Fig. 13), the SYSU and Gou and Soja (2024) products perform particularly well, with median values of 0.71 and 0.68, respectively. By comparison, WGHM performs poorly, with a median NSE of only 0.14, and GLWS shows little improvement over WGHM, with the same median NSE of 0.14. The CC shows a consistent pattern (Table 1 and Fig. 14). SYSU and Gou and Soja (2024) achieve median CC values of 0.90 and 0.88, respectively, followed by GLWS (0.68), whereas WGHM is lowest (0.65). Low NSE and CC values are concentrated in arid regions, such as northern Africa and Mongolia, likely reflecting weak hydrological signals and the limited sensitivity of both GRACE/GFO and WGHM in such environments (Gou and Soja, 2024; Xiong et al., 2025a).

To further examine overall agreement, we concatenated the basin-average time series from the 288 basins into column vectors and constructed scatter-density plots for each product, using GRACE/GFO as the reference on the x axis and the corresponding basin-average estimates derived from each product on the y axis. We then computed the R² for each product. As shown in Fig. 15, SYSU achieves the highest R² (0.92), followed by Gou and Soja (2024) (0.83), GLWS (0.45), and WGHM (0.26). Together with the results in Figs. 12–14, these comparisons indicate that all three downscaled products improve WGHM-simulated TWSA estimates to varying degrees, yielding basin-wise signals that more closely match GRACE/GFO observations, with SYSU showing the strong overall performance.

Figure 15 Comparison of three downscaled products and the WGHM output with GRACE/GFO observations at the basin scale. The color density of the scatter points represents the frequency of data point distribution. The R² characterizes the goodness of fit between each product and the observations.

In addition, another widely used assimilation-based product is the GLDAS CLSM DA v2.2 daily dataset at 0.25° resolution, which assimilates CSR mascon observations into the Catchment Land Surface Model (CLSM) using an EnKF method (Li et al., 2019a). Comparison with this product, the results show that SYSU also agrees more closely with GRACE/GFO at the basin scale. The detailed comparison is presented in the Supplement (Fig. S1 and Table S1).

5 Discussion

5.1 Uncertainty of SYSU TWSA product

The uncertainty of SYSU TWSA arises from several sources, including GRACE/GFO measurement errors, the choice of WGHM-derived spatial basis functions, ICA decomposition, and the definition of supplemental glacier and lake mascon groups. A complete propagation of all these uncertainties is difficult because WGHM does not provide explicit uncertainty estimates for its spatial structures, and the uncertainties associated with ICA mode selection and mascon group definitions are not directly available. Therefore, we provide a posterior formal uncertainty estimate from the joint-inversion adjustment under fixed spatial basis functions. This estimate reflects the uncertainty of the retrieved temporal coefficients and their propagation to the reconstructed TWSA fields, but it should not be interpreted as a complete total uncertainty budget.

The resulting uncertainty map shows clear spatial heterogeneity (Fig. 16). Most land areas have relatively low formal uncertainty, generally below 1–2 cm EWH. Higher uncertainties occur in regions with stronger hydrological variability or more complex storage processes, including parts of North America, northern India, eastern China, central and eastern Africa, Southeast Asia, and Australia. These areas require more cautious interpretation, especially where groundwater depletion, surface-water changes, or model structural limitations may affect the spatial allocation of reconstructed signals.

Figure 16 Spatial distribution of posterior formal uncertainty of SYSU TWSA. The posterior formal uncertainty was derived from the covariance matrix of the estimated temporal coefficients in the joint-inversion adjustment and propagated to the downscaled TWSA fields under fixed spatial basis functions. Units are cm equivalent water height.

5.2 Sensitivity to WGHM-derived spatial priors

The sensitivity experiment using WGHM outputs neglecting direct human impacts (WGHM-NHI) and WGHM outputs including direct human impacts (WGHM-HI) further clarifies the role of WGHM-derived spatial priors in the joint-inversion framework. As shown in Fig. S2, when WGHM outputs neglect direct human impacts, several human-induced groundwater depletion patterns are weakened or missing, especially in groundwater-dominated and heavily managed regions. By contrast, WGHM-HI better represents these depletion signals, and the corresponding joint-inversion reconstruction better preserves negative trends in regions such as Northwestern India and the North China Plain. This indicates that the regional allocation of groundwater-related signals in SYSU remains sensitive to the quality of the WGHM-derived spatial priors.

This dependence on WGHM spatial priors also explains why supplemental basis functions are needed in regions where WGHM does not explicitly represent key storage components. As shown in Fig. 11 in the main text, WGHM fails to capture glacier-related mass loss in regions such as Alaska and the Tibetan Plateau, whereas SYSU better represents these signals by incorporating dedicated glacier mascon basis functions. Similarly, the inclusion of lake mascon groups improves the representation of selected large-lake storage variations, such as those in the North American Great Lakes and Lake Victoria. These results indicate that supplemental glacier and lake basis functions can partly alleviate the limitations of WGHM-derived spatial priors in regions where the corresponding processes are missing or poorly represented.

Nevertheless, the comparison in Fig. S3 shows that the joint-inversion framework still substantially improves basin-scale consistency with GRACE/GFO, even when WGHM-NHI is used as the spatial prior. Relative to the original WGHM-NHI outputs, the joint-inversion reconstruction shows lower RMSE and higher NSE and CC across most basins. This demonstrates that the framework effectively combines the high-resolution spatial information from WGHM with the reliable large-scale temporal variability constrained by GRACE/GFO. The Figs. 11, S2 and S3 suggest that WGHM spatial priors mainly affect the regional allocation of reconstructed signals, whereas GRACE/GFO constraints improve their large-scale temporal consistency. Therefore, when the model-derived spatial priors more adequately represent the hydrological processes, or when missing processes are supplemented by dedicated basis functions, the downscaled product is expected to provide a closer approximation to the true regional distribution of TWSA.

6 Data availability

The SYSU TWSA v1.0 dataset is openly available through https://doi.org/10.11888/Terre.tpdc.303322 (Xiong et al., 2026). The GRACE/GFO Level-2 data (downloaded from https://icgem.gfz-potsdam.de/sl/temporal, last access: 1 July 2026) and mascon solution (available at http://grace.jpl.nasa.gov/data/get-data/jpl_global_mascons, last access: 1 July 2026) are given by ICGEM and JPL, respectively (https://doi.org/10.5067/GFL20-MC062, GRACE-FO, 2023; https://doi.org/10.5067/TEMSC-3JC634, Wiese et al., 2024). Monthly WGHM outputs can be accessed via https://gude.uni-frankfurt.de/home (last access: 1 July 2026). GLWS data can be publicly obtained from https://doi.org/10.1594/PANGAEA.954742 (Gerdener et al., 2023b). The downscaling product of Gou and Soja (2024) is available at https://doi.org/10.3929/ethz-b-000648738 (Gou and Soja, 2023, 2024). The GLDAS CLSM DA 2.2 product is available at https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_CLSM025_DA1_D.2.2/ (last access: 1 July 2026) (Li et al., 2019a). The ERA5-Land products are available at https://doi.org/10.24381/cds.68d2bb30 (Copernicus Climate Change Service, 2019). The level-3 and level-4 basins boundary data are available at https://www.hydrosheds.org/products/hydrobasins (last access: 1 July 2026) (Lehner and Grill, 2013). GRDC (2020) datasets are available at https://mrb.grdc.bafg.de/ (last access: 1 July 2026) (GRDC, 2020). The in situ groundwater well observations provided by the IGRAC can be downloaded at https://doi.org/10.58154/6Z0Y-DA34 (IGRAC, 2024).

7 Conclusion

We generate SYSU TWSA, a global monthly terrestrial water storage anomaly dataset at 0.5° spatial resolution spanning April 2002 to December 2022, developed to address two persistent limitations of current global gravimetry-derived downscaled products. First, publicly available global high-resolution TWSA datasets remain scarce. Second, many existing downscaling approaches rely heavily on hydrological models and therefore tend to perform poorly in regions where key GRACE/GFO observable mass variations are not adequately simulated or explicitly represented, particularly mountain glaciers and large or rapidly changing lakes. The SYSU TWSA is generated using a joint-inversion spatial downscaling framework that integrates (i) the large-scale observational constraint from GRACE/GFO, (ii) high-resolution spatial patterns from the WaterGAP Global Hydrological Model (WGHM), and (iii) self-defined mascon groups to explicitly represent mass variations associated with global mountain glaciers and selected large or rapidly changing lakes that are observable by GRACE/GFO but are not adequately represented in WGHM. By design, the framework constrains the temporal evolution of each spatial basis function using GRACE/GFO observations, while retaining high-resolution spatial detail in the final reconstruction, thereby reducing reliance on model-simulated temporal variation.

Comprehensive checks and independent comparisons demonstrate that SYSU TWSA preserves the large-scale fidelity of GRACE/GFO and provides finer spatial details. After matching effective spatial resolution, SYSU TWSA shows generally good basin-wise consistency with GRACE/GFO across basin-size classes (R²>0.85). In basins below the GRACE/GFO effective-resolution limit, SYSU-derived TWF estimates exhibit improved consistency with water-balance closure relative to raw GRACE/GFO (NSE increase of 17.1 % across 1200 level-4 basins). Independent comparisons with 28 248 in situ groundwater wells further indicate that the enhanced spatial detail is more consistent with local groundwater level variations, with correlations increasing at 67.7 % of sites relative to raw GRACE/GFO. Comparisons against representative assimilation-based and deep-learning downscaled products in both the spectral and spatial domains show that SYSU TWSA achieves competitive overall performance while providing improved representation of glacier- and lake-related signals.

SYSU TWSA is intended for applications that require a globally consistent, gravimetry-constrained TWSA dataset with enhanced spatial detail, including basin-scale water-budget studies, groundwater and cryosphere-related analyses, and intercomparison of global water storage products. Users should note that, although the framework strengthens the representation of glacier- and lake-related signals and enhances spatial detail, uncertainties remain due to GRACE/GFO measurement noise. The enhanced spatial detail mainly comes from WGHM-derived spatial patterns and the self-defined glacier and lake mascon groups, while the temporal evolution of these spatial patterns is constrained by GRACE/GFO observations. Therefore, SYSU TWSA should be interpreted as a GRACE/GFO-constrained reconstruction informed by high-resolution spatial priors, rather than as independent GRACE/GFO observations at the 0.5° grid-cell scale. Any structural biases in the WGHM spatial patterns or in the self-defined mascon groups may therefore propagate into the product. In addition, water balance-based evaluations may be affected by uncertainties in precipitation, evapotranspiration, and runoff estimates and by human interventions that are not fully represented in the basin closure framework. Similarly, the groundwater-well comparison is based on correlations between TWSA and groundwater levels and therefore reflects consistency with groundwater-level variability rather than a direct validation of groundwater storage changes.

The SYSU TWSA dataset is openly available through https://doi.org/10.11888/Terre.tpdc.303322 (Xiong et al., 2026). We anticipate that SYSU TWSA will facilitate global and regional investigations of terrestrial water storage variability by providing a high-resolution product that remains consistently constrained by satellite gravimetry observations.

Appendix A: The self-defined mascon groups of mountain glaciers and lakes

This Appendix provides the spatial configuration of the self-defined mascon groups used as supplemental spatial basis functions in the joint-inversion framework (Sect. 3.2 and 3.3). The figures are included for reference and visualization of the group definitions.

Figure A1 Global distribution of mountain glacier mascon groups. A total of 14 regions are outlined with dashed lines, each consisting of subregions represented by different colors.

Figure A2 Spatial distribution of lake Mascon groups.