Assessment of pluri-annual and decadal changes in terrestrial water storage predicted by global hydrological models in comparison with the GRACE satellite gravity mission

The GRACE (Gravity Recovery And Climate Experiment) satellite gravity mission enables global monitoring of the mass transport within the Earth's system, leading to unprecedented advances in our understanding of the global water cycle in a changing climate. This study focuses on the quantification of changes in terrestrial water storage with respect to the temporal average based on an ensemble of GRACE solutions and two global hydrological models. Significant changes in terrestrial water storage are detected at pluri-annual and decadal timescales in GRACE satellite gravity data that are generally underestimated by global hydrological models though consistent with precipitation. The largest differences (more than 20 cm in equivalent water height) are observed in South America (Amazon, São Francisco and Paraná River basins) and tropical Africa (Congo, Zambezi and Okavango River basins). Smaller but significant (a few centimetres) differences are observed worldwide. While the origin of such differences is unknown, part of it is likely to be climate-related and at least partially due to inaccurate predictions of hydrological models. Pluri-annual to decadal changes in the terrestrial water cycle may indeed be overlooked in global hydrological models due to inaccurate meteorological forcing (e.g. precipitation), unresolved groundwater processes, anthropogenic influences, changing vegetation cover and limited calibration/validation datasets. Significant differences between GRACE satellite measurements and hydrological model predictions have been identified, quantified and characterised in the present study. Efforts must be made to better understand the gap between methods at both pluri-annual and decadal timescales, which challenges the use of global hydrological models for the prediction of the evolution of water resources in changing climate conditions.

1 Introduction

The GRACE (Gravity Recovery And Climate Experiment; Tapley et al., 2004) and GRACE Follow-On (GRACE-FO; Landerer et al., 2020) missions provide spatio-temporal observations of the gravity field spanning 2 decades, sensitive to the redistribution of masses from the deep Earth's interior to the top of the atmosphere (e.g. Chen et al., 2022). The GRACE and GRACE-FO satellite observations have been widely used to estimate changes in terrestrial water storage (TWS), expressed in equivalent water heights, representing changes in surface density (i.e. changes in mass per unit area) modelled as a layer of water of variable thickness in space and time (e.g. Wahr et al., 1998). Changes in TWS range from a few millimetres to a few tens of centimetres from arid (e.g. deserts) to humid (e.g. tropical rain forests) regions of the world and are dominated first by seasonal changes and then by long-term changes including both linear trends and interannual variability (e.g. Humphrey et al., 2016). Locally (mostly along the Amazon River), seasonal TWS variations can reach up to 1 or 2 m. Decadal trends in TWS have been attributed to climate variability (e.g. change in precipitation), direct human impacts (e.g. irrigation) and the combination of both effects (Rodell et al., 2018). Significant groundwater depletion has for example been observed in the Central Valley (California), in response to two extreme and prolonged droughts intensified by groundwater pumping for agriculture, wetland management and domestic use (e.g. Scanlon et al., 2012a; Ojha et al., 2018).

Trends in TWS are often temporary due to climate variability (e.g. Alam et al., 2021) and changes in water consumption policies (e.g. Bhanja et al., 2017). Significant interannual TWS variations detected in large river basins have been attributed to a combination of eight major climate modes, including the El Niño–Southern Oscillation (ENSO), Pacific Decadal Oscillation (PDO), North Atlantic Oscillation (NAO), Atlantic Multidecadal Oscillation and Southern Annular Mode (e.g. Pfeffer et al., 2022). Successive droughts and floods events have been associated with a succession of positive (El Niño) and negative (La Niña) phases of ENSO in various regions of the world, such as Australia, southern Africa or parts of the Amazon River basin (e.g. Ni et al., 2018; Anyah et al., 2018; Xie et al., 2019). Drought (e.g. Thomas et al., 2017) and flood potential (e.g. Sun et al., 2017) indices using GRACE observations have been developed to monitor the impact of extreme events on freshwater resources, taking into account all climatic and anthropogenic mechanisms and all water reservoirs from the surface to deep aquifers.

If the spatial and temporal variability of TWS is generally well captured, global hydrological models and land surface models tend to underestimate the amplitude of seasonal signals (e.g. Döll et al., 2014a, b) and decadal trends (e.g. Scanlon et al., 2018) when compared to GRACE observations. The differences in TWS between satellite gravity observations and model predictions have been shown to depend on the choice of models and river basin considered (e.g. Döll et al., 2014b; Wada et al., 2014; Scanlon et al., 2018, 2019; Decharme et al., 2019; Yang et al., 2020; Felfelani et al., 2017). Seasonal changes in TWS are often underestimated by hydrological and land surface models in tropical, arid and semi-arid basins and overestimated at higher latitudes in the Northern Hemisphere, likely due to insufficient surface and ground water storage estimates in tropical basins and to a misrepresentation of evapotranspiration and snow physics at higher latitudes (Scanlon et al., 2019). Some models lead to better performance in heavily managed river basins and, on the contrary, to erroneous trends and seasonal cycles in regions where the natural variability is dominant (e.g. Wada et al., 2014; Scanlon et al., 2019; Felfelani et al., 2017). The performance of models also varies during recharge and discharge periods, suggesting that some processes (e.g. reservoir operation) may be adequately captured by a model, while other processes (e.g. groundwater dynamics) may be overlooked (Felfelani et al., 2017). The reasons for discrepancies between models and satellite gravity observations remain largely unknown, though improvements in the parameterisation of global hydrological and land surface models are often recommended to reliably predict spatial and temporal changes in TWS, especially regarding aquifers (e.g. Decharme et al., 2019; Scanlon et al., 2019; Felfelani et al., 2017).

This study focuses on the comparison of two global hydrological models, ISBA-CTRIP (Interaction Soil Biosphere Atmosphere – CNRM (Centre National de Recherches Météorologiques) version of Total Runoff Integrating Pathways; Decharme et al., 2019) and WGHM (WaterGap Global Hydrological Model; Müller Schmied et al., 2021), against GRACE-based TWS observations at interannual and decadal timescales. These two models have been chosen because they provide a very precise representation of hydrological processes in natural (ISBA-CTRIP) and anthropised (WGHM) environments. Besides, both models have been widely used by the scientific community. In particular, ISBA-CTRIP is contributing to the Coupled Model Intercomparison Project CMIP6 (Voldoire et al., 2019) and WGHM to the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP; Herbert and Döll, 2019). We focus here on global hydrological models, rather than land surface models using a much simpler representation of hydrological processes across continental areas. Land surface models, such as Global Land Data Assimilation System (GLDAS, Rodell et al., 2004) NOAH v3.3, generally do not take into account lateral fluxes, surface water or groundwater compartments. Such shortcomings result in less accurate estimates of TWS changes, as shown in the Supplement Sect. S1. However, land surface models, such as GLDAS NOAH, usually provide TWS estimates in near real time or at least with shorter delays than global hydrological models, making these tools essential for many hydrological applications.

While the seasonal variations in TWS have been extensively studied (e.g. Döll et al., 2014a; Wada et al., 2014; Scanlon et al., 2019; Decharme et al., 2019; Felfelani et al., 2017), little attention has been paid to longer timescales, often only estimated as linear trends (Scanlon et al., 2018; Felfelani et al., 2017). Significant non-linear variability occurs however at interannual timescales, which may lead to considerable stress on water resources and large uncertainties on climate model projections. Besides, the same model may have different performances at seasonal, interannual and decadal timescales, as different processes prevail at such different timescales (e.g. Scanlon et al., 2018, 2019; Felfelani et al., 2017). This study will therefore quantify and characterise the amplitude of TWS at interannual and decadal timescales for nine GRACE solutions (three mascon solutions and six spherical harmonic solutions) and two global hydrological models between April 2002 and December 2016.

2 Methods

2.1 Satellite gravity data

TWS changes have been estimated using the latest release of three mascon solutions from the Jet Propulsion Laboratory (JPL) (RL06 Version 02, Wiese et al., 2019), Center for Space Research (CSR) (RL06 V02; Save et al., 2016; Save, 2020) and Goddard Space Flight Center (GSFC) (RL06 V01, Loomis et al., 2019a) and six solutions based on spherical harmonic coefficients of the gravitational potential from the JPL (RL06, GRACE-FO, 2019a; Yuan, 2019), CSR (RL06, GRACE-FO, 2019b; Yuan, 2019), GeoForschungsZentrum (GFZ) (RL06, Dahle et al., 2018), Institute of Geodesy at Graz University of Technology (ITSG) (GRACE2018, Mayer-Gürr et al., 2018), Combination Service for Time-variable Gravity Fields (COST-G) (RL01, Meyer et al., 2020) and Centre National d'Etudes Spatiales - Groupe de Recherche de Géodésie Spatiale (CNES-GRGS) (RL05, Lemoine and Bourgogne, 2020). The same corrections for the geocentre (Sun et al., 2016), C₂₀ coefficients (Loomis et al., 2019b) and glacial isostatic adjustment (GIA; ICE6G-D by Peltier et al., 2018) have been applied for mascon and spherical harmonic solutions. The Stokes coefficients from the JPL, CSR, GFZ, ITSG, COST-G and CNES-GRGS solutions, with the aforementioned corrections applied, have been truncated at degree 60, converted to surface mass anomalies expressed as equivalent water height (cm) and projected onto the WGS84 ellipsoid using the locally spherical approximation (Eq. 27 in Ditmar, 2018) implemented in the l3py Python package (Akvas, 2018). Systematic errors (i.e. stripes) have been removed from spherical harmonic solutions (except for the constrained CNES-GRGS solutions) using an anisotropic filter based on the principle of diffusion (Goux et al., 2023), using Daley length scales of 200 and 300 km in the north–south and east–west directions and a shape of Matérn function close to a Gaussian (eight iterations). The diffusive filter allows for the conservation of mass within the continental domain, defined here as grid cells where at least 30 % of the altitudes from ETOPO1 Global Relief Model (NOAA National Geophysical Data Center, 2009) are above sea level. Small islands (<100 000 km²) have been excluded from the continental domain because of the limited spatial resolution of monthly GRACE products (a few hundred kilometres). By default, the GRACE-derived TWS anomaly used in this study is the average of the nine processed GRACE solutions. The uncertainty on GRACE-based TWS anomalies is estimated as the dispersion (minimum to maximum) between the nine GRACE solutions.

2.2 Global hydrological models

TWS changes have also been estimated using the ISBA-CTRIP global land surface modelling system (Decharme et al., 2019) and the 2.2d version (Müller Schmied et al., 2021) of the WGHM including glaciers.

ISBA solves the water and energy balance in the soil, canopy, snow and surface water bodies, and CTRIP simulates discharges through the global river network, as well as the dynamic of both the seasonal floodplains and the unconfined aquifers. ISBA and CTRIP are coupled through the land surface interface SURFEX, allowing complex interactions (e.g. floodplain free-water evaporation and upwards capillarity fluxes between groundwaters and superficial soils) between the atmosphere, land surface, soil and aquifer. ISBA-CTRIP is forced at a 3 h time step with the ERA-Interim atmospheric reanalysis (Dee et al., 2011) for air temperature and humidity, wind speed, surface pressure and total radiative fluxes and with the gauge-based Global Precipitation Climatology Center (GPCC) Full Data Product V6 (Schneider et al., 2014) for precipitation.

WGHM2.2d simulates changes in water flows and storage using a vertical mass balance for the canopy, snow and soil and a lateral mass balance for the surface water bodies and groundwater (Müller Schmied et al., 2021). WGHM is coupled with a global water use model, taking into account water impoundment in artificial reservoirs and regulated lakes and water withdrawals for irrigation, livestock, domestic use, manufacturing and thermal power (Müller Schmied et al., 2021). Anthropogenic water withdrawals/impoundments are assumed to only impact surface waters and groundwaters (Müller Schmied et al., 2021). In addition, water storage changes in continental glaciers have been simulated with the Global Glacier Model (Marzeion et al., 2012) and added as an input to WGHM (Cáceres et al., 2022). The WGHM uses meteorological input data from WFDEI (Weedon et al., 2014) also based on the ERA-Interim atmospheric reanalysis for air temperature and solar radiation and GPCC for precipitation. Two model variants are available using different irrigation efficiencies (optimal and 70 % of optimal) (Döll et al., 2014b). Both being equally plausible given the limited datasets available to characterise groundwater abstractions for irrigation, we averaged the two variants in the present study.

2.3 Lake data

Lake water storage anomalies were then added to the predicted TWS anomalies from ISBA-CTRIP and WGHM. Indeed, although WGHM2.2d includes artificial and natural lakes in its framework, large differences were observed between the observed and predicted TWS anomalies around large lakes (e.g. American and African Great Lakes, Caspian Sea, Volta Lake), which were greatly reduced with the application of a lake correction (Appendix A).

Changes in lake volume were estimated for 100 lakes during the whole GRACE period from the hydroweb database (https://hydroweb.theia-land.fr/, last access: 18 October 2023), based on a combination of lake level measurements from satellite altimetry and lake area measurements from satellite imagery (e.g. Crétaux et al., 2016). Then lake volume changes are converted into equivalent water heights (m) over a regular $\mathrm{1}\times \mathrm{1}{}^{\circ }$ grid, using the GLWD (Global Lakes and Wetlands Database) shapes for lakes larger than 5000 km² .

2.4 Precipitation data

Precipitation is estimated using two distinct products, the Global Precipitation Climatology Center (GPCC) Full Data Product V6 (Schneider et al., 2014) and the IMERG (Integrated Multi-satellitE Retrievals for GPM) data product (Huffman et al., 2019). GPCC is a gauge-based product. IMERG is based on the TRMM (Tropical Rainfall Measuring Mission: 2000–2015) and GPM (Global Precipitation Measurement: 2014–present) satellite data.

2.5 Data processing

The period of common availability for all datasets spans April 2002 (first estimation of TWS changes with GRACE data) to December 2016 (latest estimation of TWS changes with WGHM data). All time series have been averaged monthly. Months with missing data are excluded from all datasets, leaving 141 valid months between April 2002 and December 2016. All datasets were interpolated to a regular $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid using the conservative algorithm from xESMF (Zhuang et al., 2020), allowing the integral of the surface mass anomalies to be preserved across the grid conversion (i.e. the water mass anomaly over a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell is equal to the area-weighted average of the mass anomalies from overlapping cells in the source grid). Because this study focuses on interannual to decadal changes in total terrestrial water storage, regions where observed mass changes are known to be dominated by other processes have been masked. These include the oceans; ice-covered regions such as Antarctica, Greenland and Arctic islands; and regions impacted by very large earthquakes (Sumatra, Tohoku, Maule), defined by Tang et al. (2020). Seasonal signals have been removed by least-squares adjustment of annual and semi-annual sinusoids. Finally, to be able to compare higher-resolution hydrology products to GRACE-based TWS anomalies, a diffusive filter with an isotropic Daley length of 250 km has been applied to all products. In the following, we refer to the fully processed time series as TWS anomalies. Residual TWS anomalies (sometimes shortened as residuals) refer to the difference between the TWS anomalies estimated with the average GRACE solution and the TWS anomalies estimated with one of the two global hydrological models considered in this study (either ISBA-CTRIP or WGHM). The amplitude of the interannual variability is expressed as the range at 95 % CL (confidence level) of fully processed TWS anomalies. The range at 95 % CL is defined as the difference between the 97.5 and 2.5 percentiles. It provides a more accurate quantification of the amplitude of the non-seasonal TWS variations than the root mean square (rms), while allowing for the removal of extreme values.

Figure 1 Comparison of TWS anomalies estimated from an ensemble of nine GRACE solutions and two global hydrological models. The amplitude of the non-seasonal TWS variability is expressed as the range at 95 % CL, calculated as the difference between the 97.5 and 2.5 percentiles of the TWS anomalies obtained in each grid cell over the entire study period. TWS predictions from global hydrological models should not be compared with GRACE data around glaciers, identified by white contours. (a) Range of TWS anomalies estimated as the average of nine GRACE solutions. (b) Dispersion of the range of TWS anomalies among nine GRACE solutions. Range of TWS anomalies estimated with ISBA-CTRIP (c) and WGHM (d). Range of residual TWS anomalies estimated as the difference between the average of nine GRACE solutions and ISBA-CTRIP (e) or WGHM (f).

3 Results

3.1 Comparison of observed and predicted TWS anomalies

TWS anomalies are globally lower in hydrological models than in GRACE solutions, leaving large residuals in GRACE satellite data (Fig. 1). The underestimation of TWS anomalies is more acute with WGHM (Fig. 1d) than with ISBA (Fig. 1c). Significant (>5 cm) residual TWS anomalies (Fig. 1e and f) are observed in South America (Amazon, Orinoco, São Francisco and Paraná River basins), Africa (Congo and Zambezi basins), Australia (northern part of the continent), Eurasia (India, North European Plain, Ural Mountains, Siberian Plateau) and North America (Colorado Plateau, Rocky Mountains). All GRACE solutions are remarkably consistent with each other, which is evidenced by small dispersion values (Fig. 1b). The amplitude of non-seasonal TWS signals is very similar in mascons and spherical harmonic solutions, which is generally larger than in global hydrological models (Supplement Figs. S2.1 and S2.2).

In most regions of the world, the differences between GRACE and global hydrological models (Fig. 1e and f) are much larger than the dispersion between the different GRACE solutions. Indeed, the residual TWS anomalies are significantly larger (5th, 50th and 95th percentiles of the rms of residual TWS anomalies at 4, 8 and 20 cm) than the uncertainty on GRACE data estimated by the dispersion among the nine solutions (5th, 50th and 95th percentiles of the standard deviation between the nine GRACE solutions at 1, 3 and 13 cm). The largest (≥5 cm) dispersion values are observed in coastal and mountainous regions or in regions with very large (≥20 cm) residuals (Fig. 1b). Larger sources of errors are indeed expected near the coast in GRACE measurements due to leakage errors, making the interpretation of residual signals difficult in islands such as Madagascar or the Indonesian archipelago. Similarly significant ice melt from glaciers occurs in mountainous regions such as the Alaska or Tibetan Plateau, which is monitored by GRACE but not simulated by global hydrological models, leaving large TWS residuals (30 cm) around glaciers. Global hydrological models should therefore not be compared with GRACE around glaciers, whose limits have been determined with the sixth version of the Randolph Glacier Inventory (RGI Consortium, 2017) identified with white contours in Fig. 1.

Figure 2 Range ratios between the average GRACE solution and the hydrological models ISBA-CTRIP (a) and WGHM (b). Determination coefficients between the average GRACE solution and the hydrological models ISBA-CTRIP (c) and WGHM (d). Regions where the coefficient of determination is negative are shown in white.

To be able to differentiate a systematic underestimation of TWS anomalies from singular differences in the spatial and temporal variability, we computed the range ratio between the average GRACE solution and each hydrological model. For most regions of the world (Fig. 2a and b), the range of TWS anomalies is larger for GRACE than for ISBA-CTRIP or WGHM, except in east Canada (Ontario, Quebec, Newfoundland), north Asia (east Siberia, Ob River, Finland/north-west Russia) and central Africa (Cameroon, Gabon, Congo). In these regions, the coefficient of determination (R²) between the GRACE and the hydrological models is typically negative (Fig. 2c and d), indicating that the variance of the residuals is larger than the variance of GRACE data. The global hydrological models ISBA-CTRIP and WGHM are therefore not able to predict the TWS variability estimated from GRACE satellite data in these regions.

The large residuals observed with ISBA-CTRIP in the north-west of South America (Fig. 1e) are due to differences in the spatial and temporal variability of observed and predicted TWS changes. The range of TWS variations is indeed larger for ISBA-CTRIP than for GRACE in this region. R² values are relatively high (0.5–0.9) in the north of the Amazon, indicating important similarities between GRACE and ISBA-CTRIP. Conversely, R² values are very low (<0.3) in the south of the Amazon, indicating significant differences between GRACE and ISBA-CTRIP.

The range of TWS anomalies is smaller for hydrological models than for GRACE over most of the study area (76 % for ISBA-CTRIP and 83 % for WGHM). TWS anomalies predicted by hydrological models are underestimated by at least 50 % over almost half of the study area (40 % for ISBA-CTRIP and 49 % for WGHM). TWS anomalies are at least 2 times smaller than GRACE for 22 % of the study area for ISBA-CTRIP and 25 % for WGHM. The largest range ratios (>5) are reached across deserts (Sahara, Arabian Peninsula, Gobi Desert) and glaciers (Alaska, Patagonia, Himalaya). Such differences are due to numerical artefacts (denominator near zero) and non-hydrological signals (ice melting) observed by GRACE. Very large range ratios (2–4) are also observed for ISBA-CTRIP across the United States (Great Plains aquifer) and the north of India because of significant anthropogenic influences in these regions, with a potential contribution of glaciers across the north of India (Blazquez, 2020). Large-range ratios (from 2 to 5) are reached in tropical and subtropical regions of the Southern Hemisphere (Africa, South America, Australia) for WGHM.

Over more than half of the study area (61 % for ISBA-CTRIP and 53 % for WGHM), global hydrological models explain a minor part (R²<0.5) of the variance of the TWS anomalies estimated with the average GRACE solution (Fig. 2c and d). By comparison with GRACE, WGHM is more performant in the Northern than Southern Hemisphere. Relatively large R² values (>0.5) are reached in the United States, central and northern Europe, west and central Siberia, eastern Asia, north of India, the Caspian Sea, and the Arabian Peninsula (Fig. 2d). Large R² values are also reached over most of South America (Fig. 2d). Lower R² values (<0.5) are reached over most of the African and Australian continents and parts of the Northern (north Canada, central Asia, east Siberia, south India) Hemisphere (Fig. 2d). By comparison (Fig. 2c), ISBA-CTRIP is more performant (R²>0.5) in the Southern Hemisphere (north, central and east Australia; southern and eastern Africa; and South America except Peru, Bolivia and Patagonia) and parts of the Northern Hemisphere (eastern United States, south Canada, central and northern Europe, south of Siberia, Caspian Sea, south of India, east China). Lower R² values (<0.5) are reached for ISBA-CTRIP in north Canada, west and central Africa, Arabian Peninsula, south and central Asia, and west Australia (Fig. 2c). Both models exhibit negative R² values in central and Sahelian Africa, as well as in Quebec and Ontario (Fig. 2c and d). For ISBA-CTRIP, negative R² coefficients are also reached in north Bolivia, Alaska, north of India and Siberia (south of Lena River). For WGHM, negative R² coefficients are reached in the central United States and south India. These metrics indicate that for some regions of the world (not necessarily the same for both models), hydrological models are able to capture a large part of the TWS variability estimated from GRACE but that, overall, significant differences exist between global hydrological models and GRACE satellite data.

Figure 3 Characteristic timescales in residual TWS anomalies calculated as the differences between the average GRACE solution and ISBA-CTRIP (a) or WGHM (b). Subannual, pluri-annual and decadal contributions have been computed with high-pass (cut-off period at 1.5 years), band-pass (cut-off periods at 1.5 and 10 years) and low-pass (cut-off period at 10 years) filters respectively. The percentage of variance explained by one contribution has been calculated as the coefficient of determination with respect to the full residual signal.

3.2 Characteristic timescales of residual TWS anomalies

The differences in TWS anomalies estimated from GRACE and global hydrological models (or residual TWS anomalies) are largely dominated by pluri-annual and decadal signals (Fig. 3). Residual TWS anomalies have been separated into subannual, pluri-annual and decadal contributions using high-pass (cut-off period at 1.5 years), band-pass (cut-off periods at 1.5 and 10 years) and low-pass (cut-off period at 10 years) filters respectively. The percentage of variance explained by each contribution has been calculated as R² values and reported in Maxwell's colour triangle (Fig. 3). Residual TWS anomalies are dominated by decadal signals over a large part of the study area (51 % with ISBA-CTRIP and 40 % with WGHM), including Alaska; west Canada; the Brazilian highlands (São Francisco and Paraná River basins); Patagonia; west (Niger and Volta River basins) and southern Africa (Okavango and Zambezi River basins); parts of west (Arabic Peninsula, Caspian Sea drainage area, Tigris–Euphrates, Dnieper, Volga and Don River basins), central (Tibetan Plateau and Tarim, Ganges and Brahmaputra River basins) and north (Yenisei and Lena River basins) Asia; and east Australia. When calculating the residuals with ISBA-CTRIP, large decadal signals are also observed across north-west America (Sierra Madre, Sierra Nevada, Great Basin, Rocky Mountains) and the north of India (Indus River basin).

Pluri-annual signals are prevalent in residual TWS anomalies across central Africa, western Australia, Siberia (Ob and Yenisei), eastern Europe, north-east America (Great Lakes) and the southwest of the Amazon basin. Subannual signals are prevalent in regions with tenuous TWS variability (i.e. Sahara, southern Africa, southwest Australia), likely pointing out the remaining level of noise in GRACE data (Fig. 1b). Regions with large (≥10 cm) residual TWS anomalies (Fig. 1e) are systematically dominated by pluri-annual to decadal contributions (Fig. 3). On the other hand, regions with very small (<2 cm) residual TWS anomalies, such as the Sahara, are dominated by subannual and decadal contributions (Fig. 3). As no significant geophysical signal is expected in such regions, this can be interpreted as the spectral content of the noise, including both a high-frequency and low-frequency component.

Figure 4 (a) Linear trends in residual TWS anomalies calculated as the difference between the average GRACE solution and ISBA-CTRIP. (b) Same as (a) with WGHM. (c) Amplitude of non-linear signals in residual TWS anomalies calculated as the difference between the average GRACE solution and ISBA-CTRIP. The amplitude is calculated as the difference between the 97.5 and 2.5 percentiles. (d) Same as (c) with WGHM. (e) Coefficient of determination calculated for non-linear signals with respect to TWS anomalies calculated as the difference between the average GRACE solution and ISBA-CTRIP. (f) Same as (e) with WGHM.

Residual TWS anomalies are dominated by pluri-annual and decadal changes in the TWS, including linear trends and non-linear signals (Fig. 4). Though significant linear trends are detected (±1 cm yr⁻¹), residual TWS anomalies are mainly due to non-linear variability in the TWS (Fig. 4). Apart from glaciers, significant trends in TWS residuals are observed in west (Niger) and southern (Okavango and Zambezi) Africa, north-east Australia, south Asia (mostly the north of India, especially when using ISBA-CTRIP), north-west America (ISBA-CTRIP only) and the central United States (mainly WGHM). Part of the residual TWS trends observed with ISBA-CTRIP in north-west America (Sect. 4.6) and south Asia (Sect. 4.7) are likely due to anthropogenic influences, including groundwater abstractions primarily used for irrigation. In other regions of the world, residual trends in TWS are likely related to climate variability (precipitation excess/deficit which may be associated with the alternance of wet/dry phases of natural climate modes such as ENSO in southern Africa and north-east Australia) or land-use changes (west Africa). In most regions of the world (72 % of the study area for ISBA-CTRIP and 83 % for WGHM), the residual variability in TWS cannot be explained by a linear trend and involves significant variability at interannual and decadal timescales (Fig. 4c to f).

4 Discussion

To better characterise and understand the nature of residuals, TWS anomalies estimated from GRACE and global hydrological models have been averaged over large regions of the world and compared to in situ and satellite precipitation. In the following, we discuss regional TWS anomalies where the largest residuals are observed around the Central Amazon Corridor, the upper São Francisco River, the Zambezi and Okavango rivers, the Congo River, the north of Australia, the Ogallala aquifer in central United States, the north of the Black Sea and the Northern Plains of India (see map in Fig. B1 in Appendix B). For each of these regions, all the solutions of the GRACE ensemble (three mascon and six spherical harmonic solutions) detect slow changes in TWS, which indicates high confidence in these observations. Larger differences occur between ISBA-CTRIP and WGHM, and both models systematically underestimate the pluri-annual and decadal changes in TWS captured by GRACE. Part of these differences may be attributed to common sources of errors in GRACE-based TWS estimates, including errors in background models (for example, the atmospheric circulation model) and post-processing choices (for example, the GIA model). However, errors in the atmospheric model (GAA from AOD1B, based on ERA5) would be associated with fast changes in TWS, while errors in the GIA model (ICE6G-D) would be characterised by linear trends over the GRACE period. Here, the largest differences between GRACE and global hydrological models occur at pluri-annual and decadal timescales and are generally well correlated with precipitation. A large part of the differences between GRACE and global hydrological models are therefore likely to be climate-related and at least partially due to inaccurate predictions of global hydrological models. Similar regional analyses have been done for the 40 largest river basins of the world with comparable results (Figs. S3.1 to S3.41).

4.1 Central Amazon Corridor

4.1.1 Study area

The Central Amazon Corridor (1^∘ N–7^∘ S and 75–50^∘ W) surrounds the Solimões/Amazon mainstream river and the downstream parts of its main tributaries, including the Japura, Jurua, Purus, Negro, Madeira, Trombetas, Tapajos and Xingu rivers. Those large rivers exhibit a monomodal flood pulse lasting several months, flooding an extensive lowland area, largely covered by forests (e.g. Junk, 1997; Melack and Coe, 2021). The extension of the flooded area varies from 100 000 to 600 000 km² in the Amazon basin (e.g. Fleishmann et al., 2022), in phase with water level variations in rivers that can reach up to 15 m annually (e.g. Birkett et al., 2002; Alsdorf et al., 2007; Frappart et al., 2012; Da Silva et al., 2012), with significant interannual variability (e.g. Fassoni-Andrade et al., 2021). Heterogeneous soils distributions, including Ferralsols, Plinthosols and Gleysols (e.g. Quesada et al., 2011), lie over unconsolidated sedimentary rocks, alluvial deposits and consolidated sedimentary rocks with relatively homogeneous hydraulic properties (e.g. Gleeson et al., 2011; Fan et al., 2013). Across the central Amazon lowlands, the groundwater table fluctuates by several metres (Pfeffer et al., 2014), corresponding to groundwater storage changes of several tens of centimetres (Frappart et al., 2019), which constitutes a large part of the TWS changes observed by GRACE (Frappart et al., 2019).

Figure 5 Comparison of TWS and precipitation anomalies averaged over the Central Amazon Corridor (box A in Fig. B1 in Appendix B). (a) Average precipitation anomalies for the GPCC (gauge-based) and IMERG (satellite-based) products. (b) Power spectral density (PSD) of average precipitation anomalies. (c) TWS anomalies average over the central Amazon for two global hydrological models (ISBA-CTRIP in blue and WGHM in black) and nine GRACE solutions (mascons in red, spherical harmonic in magenta). The solid line corresponds to the average of the sub-ensemble and the shaded area to the minimum to maximum envelope. (d) PSD of the averaged TWS anomalies shown in (c). (e) Residual TWS anomalies averaged over the Central Amazon Corridor and calculated as the difference between GRACE and ISBA-CTRIP (blue when the difference is calculated with mascons and cyan with spherical harmonics) or WGHM (black when the difference is calculated with mascons and grey with spherical harmonics).

4.1.2 Comparison of global hydrological models with GRACE

Over the central Amazon region (Fig. 5), TWS anomalies predicted by global hydrological models agree well with GRACE observations, with very large Pearson coefficients reached for both ISBA-CTRIP (R=0.90) and WGHM (R=0.86). The amplitudes of TWS anomalies predicted with ISBA-CTRIP match GRACE solutions closely, while WGHM tends to underestimate the TWS variability at interannual and decadal timescales, which is likely due to a more accurate representation of the floodplains and their interactions with the atmosphere, soil and aquifer with ISBA-CTRIP than WGHM (Fig. 5d). Interannual variability occurs in the precipitation as well (Fig 5a and b), with significant correlation with GRACE (R=0.54), ISBA (R=0.59) and WGHM (R=0.64) and a phase lag of 1 month. Despite good performances for both models (especially ISBA-CTRIP), significant residual signals remain in TWS anomalies after correction of hydrological effects, consisting mostly of an increasing trend with ISBA-CTRIP, with significant interannual variability superimposed for WGHM. The residual TWS changes corrected with WGHM are still significantly correlated with precipitation (R=0.48) with a phase lag of 4 months. No significant correlation can be found between the residual TWS anomalies calculated with ISBA and precipitation anomalies (maximum R value of 0.22 with a time lag of 14 months), though significant decadal and pluri-decadal variability can be observed in GPCC precipitation records, which may explain a residual trend in TWS (∼5 mm yr⁻¹).

Residual TWS anomalies may be due to inaccurately modelled water storage variations in any reservoir from the surface to the aquifer. The largest residual TWS variations are observed along the downstream part of the Solimões, at the confluences with the Purus and the Rio Negro, which is a region that is largely covered by floodplains (e.g. Fleishmann et al., 2022) and dominated by changes in surface water storage (Frappart et al., 2019). The long timescales associated with the residuals and increasing time lags with precipitation suggest however a significant contribution from groundwater storage fluctuations, which are insufficiently constrained in global hydrological models (e.g. Decharme et al., 2019; Scanlon et al., 2018, 2019). Large floodplains may indeed delay the water transport for several months (e.g. Prigent et al., 2020), through storage and percolation from the surface towards the aquifer (e.g. Lesack and Melack, 1995; Bonnet et al., 2008; Frappart et al., 2019). Groundwater stores excess water during wet periods and sustains rivers and floodplains during low-water periods (e.g. Lesack, 1993). Groundwater systems have also been shown to convey seasonal anomalies (for example, droughts) for several years at local (e.g. Tomasella et al., 2008) and regional (Pfeffer et al., 2014) scales. Such memory effects may be underestimated by global hydrological models, which would result in much faster variations of the TWS.

4.2 Upper São Francisco

4.2.1 Study area

The São Francisco River, located in north-east Brazil, is 3200 km long and drains an area of about 630 000 km². Hydroelectric dams located along the São Francisco provide about 70 % of north-east Brazil's electricity, including the Três Marias, Sobradinho and Luíz Gonzaga (Itaparica) reservoirs with respective volumes of 15 278, 28 669 and 3549 hm³. Significant decreases in the river flow during the 1980–2015 period have been attributed to increased groundwater withdrawals sustaining irrigated agriculture and decreasing the groundwater contributions to streamflow (i.e. baseflow) (Lucas et al., 2020). As a result of a prolonged drought lasting from 2002 to 2017 (Freitas et al., 2021), the São Francisco hydroelectric plants only provided a minor part (from 18 % to 42 % depending on the year) of the total electricity demand, which was sustained by increased fossil fuel consumption (de Jong et al., 2018). A decrease in TWS was also observed from 2012 to the end of the GRACE mission (mid-2017) across the São Francisco coincident with the observed rainfall deficit (Ndehedehe and Ferreira, 2020), allowing the impact of prolonged droughts on the water supply in a vulnerable region to be better quantified (Paredes-Trejo et al., 2021).

Figure 6 Same as Fig. 5 but for the upper São Francisco (box B in Fig. B1 in Appendix B).

4.2.2 Comparison of global hydrological models with GRACE

Over the upper São Francisco region (Fig. 6), TWS anomalies predicted with global hydrological models are well correlated with GRACE data on a year-to-year basis (R=0.79 for ISBA and R=0.81 for WGHM). The times of the minimum and maximum TWS anomalies are well picked up by satellite gravity observations and models, though the amplitude of TWS anomalies is underestimated by global hydrological models. All nine GRACE solutions exhibit interannual and decadal variability in TWS, which is absent in both global hydrological models. In particular, GRACE monitors a drop in terrestrial water storage from 2012 to 2016 (Fig. 6b), corresponding to 4 years of consecutive deficit in precipitation (Fig. 6a), which is not picked up by global hydrological models. As a consequence, residual TWS anomalies (Fig. 6e), characterised by prominent interannual and decadal signals (Fig. 6f), reach 10–20 cm in the São Francisco region. TWS anomalies predicted by hydrological models are relatively well correlated with precipitation (R=0.6 for ISBA and 0.52 for WGHM) with a time lag of 1 month, while the correlation with GRACE TWS anomalies is more marginal (R=0.39 with a time lag of 1 month). Residual TWS anomalies are also only marginally correlated with precipitation (R=0.29 for GRACE-WGHM and 0.33 for GRACE-ISBA), with a time lag of 3 months.

These results tend to show that global hydrological models reproduce the year-to-year variability of TWS anomalies across the São Francisco quite well (especially in term of occurrence of a wet/dry anomaly, as the amplitudes of the anomalies may be underestimated) but struggle to predict slower hydrological processes characterised by interannual and decadal timescales.

4.3 Zambezi–Okavango

4.3.1 Study area

The Zambezi River basin, located in south tropical Africa, drains an area of 1.4×10⁶ km² connecting Angola (18.3 %), Namibia (1.2 %), Botswana (2.8 %), Zambia (40.7 %), Zimbabwe (15.9 %), Malawi (7.7 %), Tanzania (2.0 %) and Mozambique (11.4 %) (Vörösmarty and Moore, 1991). It encompasses humid, semi-arid and arid regions dominated by seasonal rainfall patterns associated with the Intertropical Convergence Zone (ITCZ), with a wet season spanning October to April and a dry season spanning May to September (Lowmann et al., 2018). The Zambezi basin harbours very large wetland areas and lakes, whose extension considerably varies with precipitation at seasonal and interannual timescales (Hugues et al., 2020). Significant interannual variability in the precipitation and TWS has been detected over the Zambezi and Okavango regions and attributed to several climate modes, including the Pacific Decadal Oscillation, Atlantic Multidecadal Oscillation and El Niño–Southern Oscillation (Pfeffer et al., 2021).

Figure 7 Same as Fig. 5 but for the Zambezi and Okavango rivers (box C in Fig. B1 in Appendix B).

4.3.2 Comparison of global hydrological models with GRACE

Across the Zambezi and Okavango region (Fig. 7), TWS anomalies are well correlated with precipitation (R=0.62 and 0.49 with a time lag of 1 month for ISBA-CTRIP and WGHM). Positive (respectively negative) precipitation anomalies correspond to a local maximum (respectively minimum) in TWS. This year-to-year variability is consistent between GRACE and global hydrological models, as evidenced by a Pearson correlation coefficient of 0.60 between GRACE and ISBA-CTRIP and 0.63 between the GRACE and WGHM. However, the TWS anomalies estimated from GRACE exhibit a strong decadal oscillation, with a minimum in 2005/2006 and a maximum in 2011/2012 that is not picked up by hydrological models, leaving a very strong (20 cm in amplitude) decadal anomaly in the residuals TWS. Though the residual TWS anomalies are poorly correlated with the precipitation anomaly (R=0.23 and 0.25 with a phase lag of 28 and 40 months for GRACE–ISBA and GRACE–WGHM respectively), they are strongly related to the accumulated precipitation anomalies, also exhibiting a strong decadal anomaly with a minimum in 2005/2006 and a maximum in 2011/2012.

The TWS residuals can be reduced locally by up to 50 % in the Zambezi region by applying an empirical model based on climate modes, as formulated by Pfeffer et al. (2021). The main modes of variability found in the TWS residuals are the Pacific Decadal Oscillation and the Atlantic Multidecadal Oscillation.

4.4 Congo

4.4.1 Study area

The Congo basin is the second-largest river basin in the world, with a drainage area of $\sim \mathrm{3.7}\times {\mathrm{10}}^{\mathrm{6}}$ km² and an average annual discharge of ∼40 500 m³ s⁻¹ (Laraque et al., 2020). Despite its importance, the Congo River basin is scarcely studied (Alsdorf et al., 2016), though a growing interest has arisen over the past decade, substantially due to advances in satellite hydrology (e.g. Papa et al., 2022; Paris et al., 2022; Schumann et al., 2022). With an average rainfall around 1500 mm⁻¹, the Congo basin benefits from a humid tropical climate with a complex seasonal migration of rainfall across the basin with a first maximum in November–December and a second peak in April–May (Alsdorf et al., 2016), leading to a bimodal river discharge (Kitambo et al., 2022). The “Cuvette centrale” is a topographic depression located at the centre of the basin, harbouring wetlands covered by rainforests permanently or periodically flooded (Becker et al., 2018). The Congo floodplain hydrodynamics are disconnected from the main river, with much less variability observed throughout the year (Alsdorf et al., 2016). The Congo River basin hosts a large complex fractured sedimentary aquifer, with relatively low storage but high recharge rates (Scanlon et al., 2022). Very little is known about the groundwater storage variability, though comparisons of satellite estimations of the surface water storage with the total terrestrial water storage changes from GRACE suggest that most (∼90 % at annual timescales) of the variability in water storage occurs under the surface (Becker et al., 2018).

Figure 8 Same as Fig. 5 but for the Congo River (box D in Fig. B1 in Appendix B).

4.4.2 Comparison of global hydrological models with GRACE

Non-seasonal TWS anomalies are very different over the Congo basin depending on the method of estimation considered (Fig. 8). All nine GRACE solutions are consistent with each other but differ from both global hydrological models that also exhibit large discrepancies with each other (Fig. 8). The correlations of TWS anomalies with precipitation are also marginal (maximum correlation of 0.5 with WGHM). All nine GRACE solutions exhibit a 6-year cycle, in phase with accumulated precipitation with local minima in 2006 and 2012 and local maxima in 2003, 2009 and 2015 (Fig. 8). Slow changes in TWS observed with GRACE are not predicted by hydrological models, leaving large residuals in TWS characterised by a ∼6-year cycle (Fig. 8).

Significant power is found in multi-decadal precipitation time series at similar periods, ranging from 5 to 8 years (Laraque et al., 2020), as well as in discharge times series at 7.5 and 13.5 years (Labat et al., 2005). The variability of the TWS cannot be explained by major climate modes over the Congo River basin, except for the PDO, which may slightly influence the TWS variability in the north of the Congo River (Pfeffer et al., 2022). The variability in river discharge has been found to be temporarily consistent with NAO at 7.5 years (from the 1970s to the 1990s) and 35 years (from the 1940s to the 1990s) (Labat et al., 2005). Part of the inaccuracies in global hydrological models may be due to (i) the scarcity of in situ data available to constrain precipitation (Fig. 2 in Laraque et al., 2020); (ii) errors in runoff and evapotranspiration fluxes; or (iii) unresolved underground processes, including preferential flow along faults (Fig. 1 in Garzanti et al., 2019).

4.5 North Australia

4.5.1 Study area

The climate of north Australia is characterised by a wet season lasting from November to April, subject to intense thunderstorms and cyclones, with virtually no precipitation during the remainder of the year (Smith et al., 2008). Annual streamflow is highly dominated by monsoon rainfall, with dry-season flows fed by groundwater discharge that may stop for several months for a large number of rivers (Petheram et al., 2008; Smerdon et al., 2012). Groundwater plays an essential role in north Australia as it sustains rivers and vegetation, through baseflow and water uptake for plant transpiration (Lamontagne et al., 2005; O'Grady et al., 2006). Significant interannual variability, principally related to ENSO in the north of the continent, has been observed in rainfall (Cai et al., 2001; Sharmila and Hendon, 2020), river discharge (Chiew et al., 1998; Ward et al., 2010) and terrestrial water storage (Xie et al., 2019). During the GRACE era, Australia encountered a prolonged drought from 2002 to 2009, sometimes referred to as the “millennium drought” or “big dry”, immediately followed by intensely wet conditions in 2010–2011 (the “big wet” associated with La Nina) and a sustained drought, leading to another dry El Niño event in 2015 (Fig. 3 in Xie et al., 2019, and Fig. 9 in the present paper). Three major climate modes (ENSO, IOD and SAM) are necessary to explain the water storage variability across Australia, but the northern part of the country is dominated by ENSO (Xie et al., 2019).

Figure 9 Same as Fig. 5 but for north Australia (box E in Fig. B1 in Appendix B).

4.5.2 Comparison of global hydrological models with GRACE

Across north Australia (Fig. 9), TWS anomalies predicted by global hydrological models are well correlated with precipitation (R=0.73 and 0.67 with a phase lag of 1 month for ISBA and WGHM) and TWS anomalies estimated with GRACE (R=0.76 and 0.71 with ISBA and WGHM respectively). The amplitude of extreme events (for example La Niña in 2011) from ISBA matches GRACE estimates, while WGHM tends to underestimate the response of TWS to both dry (2005) and wet (2011) events (Fig. 9). The main difference between TWS estimations from global hydrological models and GRACE solutions is the pace at which TWS returns to average conditions after a wet/dry event (Fig. 9). For example, after the flooding events associated with La Niña 2011, all nine GRACE solutions estimate a slow decrease in the TWS returning to average conditions in about two years (Fig. 9). On the other hand, both global hydrological models predict a sharp decrease in the TWS returning to average conditions in about 6 months (Fig. 9). As a consequence, a positive TWS anomaly remains in the residuals after La Niña (Fig. 9), accounting for the differences in the rate of change of TWS.

These results are consistent with the findings of Yang et al. (2020), who found that except for the CLM-4.5 model, hydrological models underestimated the GRACE-derived TWS trends across Australia, due to inaccurately modelled contributions from soil moisture and groundwater storage. Similarly, TWS anomalies from GRACE were found to be a better link between vegetation change and climate variability than precipitation (Xie et al., 2019) because they convey more information about water availability in the soils and aquifers, especially when associated with SMOS (Soil Moisture and Salinity) measurements (Tian et al., 2019).

4.6 Central United States: Ogallala aquifer

4.6.1 Study area

The Ogallala, or High Plains, aquifer covers a surface area of about 450 000 km² across eight states in the central United States, including parts of Colorado, Kansas, Nebraska, New Mexico, Oklahoma, South Dakota, Texas and Wyoming. The Ogallala aquifer region supports about 20 % of the wheat, corn and cotton production in the United States (Houston et al., 2013). Groundwater abstractions for irrigation began in Texas in the 1930s (Luckey et al., 1981) and exceeded recharge over much of the central and southern parts of aquifer in the 1950s (Luckey and Becker, 1999), resulting in substantial decline of the groundwater table in the Southern and Central High Plains, while the Northern High Plains stayed in balance or replenished (Haacker et al., 2016). At current depletion rates, a large part of irrigation (about 30 %) may not be supported in the coming decades (Scanlon et al., 2012b; Haacker et al., 2016; Steward and Allen, 2016; Deines et al., 2020).

Figure 10 Same as Fig. 5 but for the central United States – Ogallala aquifer region (box F in Fig. B1 in Appendix B).

4.6.2 Comparison of global hydrological models with GRACE

In the Ogallala aquifer region, all GRACE solutions exhibit a series of upwards and downwards trends in TWS with a regular increase from mid-2006 to mid-2011 and a sharp decrease in TWS from mid-2011 to 2013, followed by another increase in TWS from early 2013 to 2016 (Fig. 10). This pattern is linked with precipitation anomalies that were mainly in excess over 2006–2011 and in deficit over 2011/2013 and oscillated around average values over 2013–2016, with a remarkably rainy year in 2014 (Fig. 10). This succession of opposite trends is not predicted by global hydrological models (Fig. 10). WGHM does predict a sharp decrease in TWS from mid-2011 to 2013 but fails to predict the increase in TWS during 2006–2011 in spite of abundant precipitation (Fig. 10).

Such differences might be explained by an overestimation of water abstractions by WGHM, which would result in almost constant TWS changes, while precipitation, and subsequent aquifer recharge, is increasing. This assumption is supported by the work of Rateb et al. (2020), showing that global hydrological models such as WGHM or PCR-GLOBWB tend to overestimate groundwater depletion due to human intervention in the region. Good agreement is found between GRACE and in situ observations of the groundwater table, though large uncertainties affect (i) the decomposition of the GRACE-based TWS anomalies into individual water reservoirs (Brookfield et al., 2018) and (ii) the estimation of hydraulic parameters (i.e. conductivity and specific yield), allowing for the conversion of groundwater level variations to groundwater storage variations (Seyoum and Milewski, 2016). For the Ogallala aquifer region, GRACE data may help to characterise insufficiently well constrained parameters of WGHM, such as hydraulic parameters (i.e. conductivity, specific yield), or parameters of the water use model, such as irrigation efficiencies. In its current stage, the ISBA-CTRIP model is not adapted to estimate TWS changes in heavily managed regions because it does not take irrigation into account.

4.7 North of India

4.7.1 Study area

The north of India hosts the Indus, Ganges and Brahmaputra River basins, with an average annual rainfall of 545, 1088 and 2323 mm yr⁻¹ respectively (e.g. Bhanja et al., 2016). The average population density ranges from 26–250 persons per km² in the north-west of India to over 1000 persons per km² in the north-east of India (Dangar et al., 2021). India is the largest groundwater user in the world, with an annual withdrawal of 230 km³ for irrigation, used essentially for rice, wheat, sugarcane, cotton and maize cultures (Mishra et al., 2018; Xie et al., 2019). High abstraction rates largely exceeding precipitation rates have been reported in north-west India, in particular in the Punjab region, leading to an aquifer depletion rate of about 1 m yr⁻¹ (Mishra et al., 2018; Dangar et al., 2021). The Northern Plains of India are bordered by the Southern Tibetan Plateau, whose glaciers have been undergoing significant ice thinning due to increased temperatures (e.g. Hugonnet et al., 2021). Both contributions from land hydrology and glaciers may therefore influence GRACE-based TWS estimates in this region.

Figure 11 Same as Fig. 5 but for the Northern Plains of India (box G in Fig. B1 in Appendix B).

4.7.2 Comparison of global hydrological models with GRACE

Because WGHM takes irrigation into account, predicted TWS anomalies closely match GRACE observations (R=0.96), leaving residuals of about ±2.5 cm (Fig. 11), which is about 4 to 6 times less than across the central Amazon (Fig. 5) or Zambezi (Fig. 7) regions. As expected in strongly anthropised regions, ISBA-CTRIP fails to recover the TWS changes estimated with GRACE, characterised by a clear decreasing trend ($-\mathrm{7.71}\pm \mathrm{0.71}$ mm yr⁻¹) over 2002–2016 (Fig. 11), clearly due to groundwater abstractions for irrigation.

Besides, the superposition of several sources of mass redistributions (i.e. land hydrology and glaciers) may generate ambiguities in the interpretation of GRACE-based TWS estimates in the north of India (Blazquez, 2020). Groundwater abstractions were however found to be the dominant driver of water mass losses across northern India (e.g. Xiang et al., 2016). Numerous studies have reported a good agreement between in situ groundwater level measurements and GRACE TWS measurements in the north of India (e.g. Bhanja et al., 2016; Dangar et al., 2021). Detailed studies indicated that better model performances could be gained by adjustment of several parameters (water percolation rate, crop water stress, irrigation efficiency, soil evaporation compensation and groundwater recession) against GRACE data (Xie et al., 2019). Such information is critical to ensure the reliability of hydrological models across several regions. For example, the ISBA-CTRIP model exhibits better performances than WGHM when compared to GRACE across the southern Peninsular Plateau of India (Fig. 1) because of an overestimation of groundwater abstractions in WGHM, leading to spurious decreasing trends not observed by satellite gravity measurements (Fig. S4.1). An increase in TWS and in replenishment of groundwater resources has indeed been reported in south India from the analysis of GRACE and well data (e.g. Asoka et al., 2017; Bhanja et al., 2017).

4.8 North of the Black Sea

4.8.1 Study area

The Black Sea catchment hosts a population of 160 million people in 23 countries drained by major rivers including the Danube, Dniester, Dnieper, Don, Kuban, Sakarya and Kizirmak. The annual precipitation varies from less than 190 mm yr⁻¹ at the north-east of the catchment (Russia) to more than 3000 mm yr⁻¹ at the west (south Austria, Slovenia, Croatia) (Rouholahnejad et al., 2014, 2017). The annual average temperature varies from 2 to 7 ^∘C at the north of the catchment (East European Plain at the border of Ukraine, Belarus and Russia), with a local minimum ($<-\mathrm{3}$ ^∘C) in the Krasnodar region (southwest Russia), to over 15 ^∘C at the south of the catchment (north of Turkey) (Rouholahnejad et al., 2014, 2017). Land use in the Black Sea catchment is dominated by agriculture (Rouholahnejad et al., 2014, 2017).

Figure 12 Same as Fig. 5 but for the north of the Black Sea (box H in Fig. B1 in Appendix B).

4.8.2 Comparison of global hydrological models with GRACE

Large TWS residuals are observed in the north-east of the Black Sea catchment, in the East European Plain crossing Ukraine, Belarus and Russia (Fig. 12). Large (∼20 cm) TWS changes are observed by GRACE satellites in this region, characterised by a decreasing trend conjugated with significant interannual variability, with a peak at 6–7 years (Fig. 12). Such TWS changes are not predicted by hydrological models, leaving large (∼15 cm) TWS residuals, dominated by decadal and interannual variability (Fig. 12).

Due to rising temperatures, a generalised drop (10 %–15 %) in solid precipitation has been observed across the East European Plain, partially offset by liquid precipitation, except along the northern coast of the Black and Azov Sea (drop ∼10 %), the lower Volga River basin (drop ∼20 %) and the Dvina River basin further north (drop ∼25 %) (Kharmalov and Kireeva, 2020). A drop in summer precipitation, together with an increase in temperature, was observed at the north of the Black, Azov and Caspian Sea, generating severe drought conditions in the region (Kharmalov and Kireeva, 2020). Water scarcity has indeed become a critical concern, with increased water stress and decreased water availability, observed today and predicted to increase in the future (Rouholahnejad et al., 2014, 2017).

5 Conclusions

Over most (>75 %) continental areas, non-seasonal TWS anomalies are underestimated by the global hydrological models ISBA-CTRIP and WGHM when compared to GRACE solutions. While both hydrological models agree relatively well with GRACE observations on short timescales (i.e. typically less than 2 years), they systematically underestimate slower changes in TWS observed by GRACE satellites occurring on pluri-annual to decadal timescales. Particularly large (15–20 cm) residual TWS anomalies are observed across the north-east of South America (Orinoco, Amazon and São Francisco basins), tropical Africa (Zambezi and Congo rivers basin) and north Australia.

In such remote areas, better performances are reached with ISBA-CTRIP than WGHM, owing to the detailed representation of hydrological processes in a natural environment. However, the TWS predicted with ISBA-CTRIP still lacks amplitude at pluri-annual and decadal timescales, leaving large linear (Amazon) and non-linear (São Francisco, Zambezi, Congo, north Australia) trends in the TWS residuals.

The comparison of global hydrological models against GRACE data does not allow for the identification of the processes responsible for these discrepancies that could originate from any reservoir from the surface to deep aquifers. However, long timescales associated with the residuals, combined with increasing time lags and decreasing correlations with precipitation, suggest at least some mismodelled contributions from the groundwater cycle. Aquifers constitute the natural accumulation of runoff and precipitation, and misestimated parameters (hydraulic properties such as the conductivity or storage capacity) and flows (e.g. recharge, discharge, deep inflow, preferential flow along faults and fractures) may lead to significant errors in predicted groundwater storage changes. An overestimation of runoff and/or evapotranspiration may also lead to an excessively quick return of the water to the atmosphere and ocean. Evapotranspiration may in particular be difficult to estimate in regions with temporary surface water bodies (for example, related to the variation of the floodplain extension or to the formation of temporary rivers flowing during the wet season and dried up during the dry season).

If ISBA-CTRIP leads to TWS predictions in better agreement with GRACE than WGHM over remote areas, the situation is inverted for strongly anthropised regions such as the Northern Plains of India and the Central Valley (California, United States) or Great Plains (Ogallala, United States) aquifer regions. Unlike WGHM, ISBA-CTRIP does not account for human-induced changes in the TWS and is therefore not able to reproduce TWS changes in highly anthropised regions. However, important differences between GRACE and WGHM are still observed in some highly anthropised regions, such as the Ogallala aquifer, which may be due to locally misestimated parameters.

Large uncertainties may indeed affect the parameterisation of the water use model. For example, an overestimation of the irrigation efficiency may lead to an overestimation of evapotranspiration and underestimation of deep percolation. Errors in such parameterisation may have a strong effect on the predicted TWS changes that could eventually be more accurately estimated using GRACE to constrain unknown parameters. The calibration and evaluation of global hydrological models would therefore benefit the consideration of a broader range of datasets, including traditional discharge data but also including terrestrial water storage anomalies from GRACE satellites. For example, WGHM simulations were shown to be improved by the joint calibration against water discharge and GRACE-based TWS anomalies (Werth et al., 2009). GRACE-based observations have also been proven useful to quantify the impact of irrigation on groundwater resources in northern India and improve groundwater forecasts under different Representative Concentration Pathways (RCPs) in the region (Xie et al., 2020). The assimilation of GRACE and GRACE-FO observations into global hydrological models can also increase the model performance with various applications (see Soltani et al., 2021, for a review). Among them, GRACE data assimilation can be used to increase the accuracy (e.g. Zaitchik et al., 2008) or resolution (e.g. Kumar et al., 2016) of predicted TWS changes. GRACE data assimilation can also be used to achieve a better separation of TWS changes into the different water storage compartments (i.e. snow, canopy, surface, soil, aquifer), usually using several other remote sensing datasets to constrain the water storage changes in individual compartments (e.g. Tian et al., 2018) or taking advantage of the integrated nature of GRACE measurements to better constrain water storage compartments that are difficult to access, such as groundwater (e.g. Girotto et al., 2017; Li et al., 2019). Significant advances would be expected from the generalisation of such approaches in a dedicated framework (e.g. Condon et al., 2021; Gleeson et al., 2021).

Appendix A: Comparison of TWS anomalies from GRACE and global hydrological models over large lakes

Residual TWS anomalies (Fig. A1) are compared for ISBA-CTRIP and WGHM with and without including the lake correction from the hydroweb database based on satellite altimetry and satellite imagery measurements. The TWS residuals are reduced for both models when applying the lake correction, especially around the Caspian Sea (−30 cm), North American Great Lakes (−7 cm), African Great Lakes (−15 cm) and Volta Lake (−5 cm). A marginal increase (+2 cm) in TWS residuals can be observed for high altitude lakes of the Tibetan Plateau (e.g. Pu Moyongcuo, Yamzho Yumco, Namu Cuo, Qinghai). Slight increases in the TWS residuals (at most +1 cm) are observed in a few anthropised regions when applying the lake correction to ISBA-CTRIP, especially near the Zeya Reservoir (Russia) and the Roraima region (north Brazil). Overall, the prediction of TWS anomalies due to hydrology is improved when using the lake correction, and the residual TWS anomalies are reduced.

Figure A1 (a) Range of residual TWS anomalies calculated with ISBA-CTRIP. (b) Range of residual TWS anomalies calculated with WGHM. (c) Range of residual TWS anomalies calculated with ISBA-CTRIP without including the lake correction. (d) Range of residual TWS anomalies calculated with WGHM without including the lake correction. (d) Difference between (a) and (c) due to the lake correction. (e) Difference between (b) and (d) due to the lake correction.

Appendix B: Location of eight regions with significant residual TWS anomalies

Residual TWS anomalies are calculated as the difference between the TWS anomalies estimated from GRACE and global hydrological models. The ensemble of residual TWS anomalies counts 18 solutions, pertaining to nine GRACE solutions (three mascon and six spherical harmonic solutions) and two global hydrological models (ISBA-CTRIP and WGHM). The range of average residual TWS anomalies shown in Fig. B1a depends on the systematic biases between the TWS estimates from GRACE and global hydrological models. These differences are significant if they exceed the dispersion among the 18 solutions, calculated as the difference between the 97.5 and 2.5 percentiles of the range of residual TWS anomalies (see Fig. B1b). The significance ratio of residual TWS anomalies (Fig. B1c) has been calculated to identify where the differences between GRACE solutions and hydrological models are significant, regardless of the solution or model considered. The dispersion of residual TWS solutions (Fig. B1b) is much larger than the dispersion of GRACE-based TWS solutions (Fig. 1b), showing that the differences between the two models may have a large impact on the residuals and their significance.

To explore a large variety of scenarios, we selected eight regions with large residuals (>10 cm) and high significance ratio (>2), including the Central Amazon Corridor (region A), the upper São Francisco River (region B), the Zambezi and Okavango rivers (region C), the Congo River (region D), the north of Australia (region E), the Ogallala aquifer in central United States (region F), the north of the Black Sea (region H) and the Northern Plains of India (region G). It may be noted that the significance ratio is not extremely high across the north of India because of the differences in the predictions of ISBA-CTRIP and WGHM. Region G was included to discuss the differences between models with respect to GRACE-based TWS anomalies. Glaciers and coastal regions have been excluded from the analyses (see Sect. 3.1).

Figure B1 (a) Average range of 18 residual TWS anomalies. (b) Dispersion of the range of residual TWS anomalies. The dispersion is calculated as the difference between the 97.5 and 2.5 percentiles of the range of 18 residual TWS anomalies. (c) Significance ratio of the averaged residual TWS anomalies calculated as the average range of residual TWS anomalies (a) divided by the dispersion of the range among the 18 solutions (b).