Julia Pfeffer
email: julia.pfeffer@magellium.fr

Anny Cazenave

Alejandro Blazquez

Bertrand Decharme

Simon Munier

Anne Barnoud

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.

Introduction

The GRACE (Gravity Recovery And Climate Experiment; [START_REF] Tapley | GRACE measurements of mass variability in the Earth system[END_REF] and GRACE Follow-On (GRACE-FO; [START_REF] Landerer | Extending the global mass change data record: GRACE Follow-On instrument and science data performance[END_REF] 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. [START_REF] Chen | Applications and challenges of GRACE and GRACE follow-on satellite gravimetry[END_REF]. 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. [START_REF] Wahr | Time variability of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE[END_REF]. 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. [START_REF] Humphrey | Assessing global water storage variability from GRACE: Trends, seasonal cycle, subseasonal anomalies and extremes[END_REF]. Locally (mostly along the Amazon River), seasonal TWS variations can reach up to 1 or 2 m. Decadal trends in TWS have been Published by Copernicus Publications on behalf of the European Geosciences Union.

attributed to climate variability (e.g. change in precipitation), direct human impacts (e.g. irrigation) and the combination of both effects [START_REF] Rodell | Emerging trends in global freshwater availability[END_REF]. 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;[START_REF] Ojha | Sustained groundwater loss in California's Central Valley exacerbated by intense drought periods[END_REF].

Trends in TWS are often temporary due to climate variability (e.g. [START_REF] Alam | Post-Drought Groundwater Storage Recovery in California's Central Valley[END_REF] and changes in water consumption policies (e.g. [START_REF] Bhanja | Groundwater rejuvenation in parts of India influenced by water-policy change implementation[END_REF]. 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. [START_REF] Pfeffer | Analysis of the interannual variability in satellite gravity solutions: detection of climate modes fingerprints in water mass displacements across continents and oceans[END_REF]. 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. [START_REF] Ni | Global terrestrial water storage changes and connections to ENSO events[END_REF][START_REF] Anyah | Understanding linkages between global climate indices and terrestrial water storage changes over Africa using GRACE products[END_REF][START_REF] Xie | Multi-climate mode interactions drive hydrological and vegetation responses to hydroclimatic extremes in Australia, Remote Sens[END_REF]. Drought (e.g. [START_REF] Thomas | GRACE groundwater drought index: Evaluation of California Central Valley groundwater drought[END_REF] 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. [START_REF] Scanlon | Global models underestimate large decadal declining and rising water storage trends relative to GRACE satellite data[END_REF] 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;[START_REF] Scanlon | Global models underestimate large decadal declining and rising water storage trends relative to GRACE satellite data[END_REF][START_REF] Scanlon | Tracking seasonal fluctuations in land water storage using global models and GRACE satellites[END_REF][START_REF] Decharme | Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications[END_REF][START_REF] Yang | Spatio-temporal evaluation of water storage trends from hydrological models over Australia using GRACE mascon solutions[END_REF][START_REF] Felfelani | Natural and human-induced terrestrial water storage change: A global analysis using hydrological models and GRACE[END_REF]. 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 [START_REF] Scanlon | Tracking seasonal fluctuations in land water storage using global models and GRACE satellites[END_REF]. 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;[START_REF] Scanlon | Tracking seasonal fluctuations in land water storage using global models and GRACE satellites[END_REF][START_REF] Felfelani | Natural and human-induced terrestrial water storage change: A global analysis using hydrological models and GRACE[END_REF]. 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 [START_REF] Felfelani | Natural and human-induced terrestrial water storage change: A global analysis using hydrological models and GRACE[END_REF]. The reasons for discrepancies between models and satellite gravity ob-servations 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. [START_REF] Decharme | Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications[END_REF][START_REF] Scanlon | Tracking seasonal fluctuations in land water storage using global models and GRACE satellites[END_REF][START_REF] Felfelani | Natural and human-induced terrestrial water storage change: A global analysis using hydrological models and GRACE[END_REF].

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; [START_REF] Decharme | Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications[END_REF] 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 [START_REF] Voldoire | Evaluation of CMIP6 DECK experiments with CNRM-CM6-1[END_REF] and WGHM to the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP; [START_REF] Herbert | Global assessment of current and future groundwater stress with a focus on transboundary aquifers[END_REF]. 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, [START_REF] Rodell | The global land data assimilation system[END_REF]) 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;[START_REF] Scanlon | Tracking seasonal fluctuations in land water storage using global models and GRACE satellites[END_REF][START_REF] Decharme | Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications[END_REF][START_REF] Felfelani | Natural and human-induced terrestrial water storage change: A global analysis using hydrological models and GRACE[END_REF], little attention has been paid to longer timescales, often only estimated as linear trends [START_REF] Scanlon | Global models underestimate large decadal declining and rising water storage trends relative to GRACE satellite data[END_REF][START_REF] Felfelani | Natural and human-induced terrestrial water storage change: A global analysis using hydrological models and GRACE[END_REF]. 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. [START_REF] Scanlon | Global models underestimate large decadal declining and rising water storage trends relative to GRACE satellite data[END_REF][START_REF] Scanlon | Tracking seasonal fluctuations in land water storage using global models and GRACE satellites[END_REF][START_REF] Felfelani | Natural and human-induced terrestrial water storage change: A global analysis using hydrological models and GRACE[END_REF]. 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.

Methods

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, [START_REF] Wiese | JPL GRACE and GRACE-FO Mascon Ocean, Ice, and Hydrology[END_REF], Center for Space Research (CSR) (RL06 V02; [START_REF] Save | High resolution CSR GRACE RL05 mascons[END_REF][START_REF] Save | CSR GRACE and GRACE-FO RL06 Mascon Solutions v02[END_REF] 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;[START_REF] Yuan | GRACE Follow-On Level-2 Gravity Field Product User Handbook[END_REF], CSR (RL06, GRACE-FO, 2019b;[START_REF] Yuan | GRACE Follow-On Level-2 Gravity Field Product User Handbook[END_REF], GeoForschungsZentrum (GFZ) (RL06, [START_REF] Dahle | GRACE Geopotential GSM Coefficients GFZ RL06, V. 6.0[END_REF], Institute of Geodesy at Graz University of Technology (ITSG) (GRACE2018, [START_REF] Mayer-Gürr | ITSG-Grace2018 -Monthly, Daily and Static Gravity Field Solutions from GRACE[END_REF], Combination Service for Timevariable Gravity Fields (COST-G) (RL01, [START_REF] Meyer | International Combination Service for Time-variable Gravity Fields[END_REF] and Centre National d'Etudes Spatiales -Groupe de Recherche de Géodésie Spatiale (CNES-GRGS) (RL05, [START_REF] Lemoine | RL05 monthly and 10-day gravity field solutions from CNES/GRGS[END_REF]. The same corrections for the geocentre [START_REF] Sun | Optimizing estimates of annual variations and trends in geocenter motion and J2 from a combination of GRACE data and geophysical models[END_REF], C 20 coefficients (Loomis et al., 2019b) and glacial isostatic adjustment (GIA; ICE6G-D by [START_REF] Peltier | Comment on the paper by Purcell et al. 2016 entitled An assessment of ICE-6G_C (VM5a) glacial isostatic adjustment model[END_REF] 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 [START_REF] Ditmar | Conversion of time-varying Stokes coefficients into mass anomalies at the Earth's surface considering the Earth's oblateness[END_REF] 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 northsouth 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 2 ) 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.

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 [START_REF] Dee | The ERA-Interim reanalysis: Configuration and performance of the data assimilation system[END_REF] 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 [START_REF] Weedon | The WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied to ERA-Interim reanalysis data[END_REF] 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.

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 measure-ments from satellite altimetry and lake area measurements from satellite imagery (e.g. [START_REF] Crétaux | Lake volume monitoring from space[END_REF]. Then lake volume changes are converted into equivalent water heights (m) over a regular 1 × 1 • grid, using the GLWD (Global Lakes and Wetlands Database) shapes for lakes larger than 5000 km 2 .

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 [START_REF] Huffman | GPM IMERG Final Precipitation L3 1 month 0[END_REF]. 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.

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 1 • × 1 • grid using the conservative algorithm from xESMF [START_REF] Zhuang | JiaweiZhuang/x-ESMF: v0.3.0 Adding ESMF[END_REF], allowing the integral of the surface mass anomalies to be preserved across the grid conversion (i.e. the water mass anomaly over a 1 • × 1 • 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 [START_REF] Tang | Seismic impact of large earthquakes on estimating global mean ocean mass change from GRACE[END_REF]. 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.

Results

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 andf) 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 andf) 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.

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 hydrologi- cal model. For most regions of the world (Fig. 2a andb), 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 2 ) between the GRACE and the hydrological models is typically negative (Fig. 2c andd), 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 2 values are relatively high (0.5-0.9) in the north of the Amazon, indicating important similarities between GRACE and ISBA-CTRIP. Conversely, R 2 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 [START_REF] Blazquez | Satellite characterization of water mass exchange between ocean and continents at interannual to decadal timescales[END_REF]. 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 2 < 0.5) of the variance of the TWS anomalies estimated with the average GRACE solution (Fig. 2c andd). By comparison with GRACE, WGHM is more performant in the Northern than Southern Hemisphere. Relatively large R 2 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 2 values are also reached over most of South America (Fig. 2d). Lower R 2 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 2 > 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 2 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 2 values in central and Sahelian Africa, as well as in Quebec and Ontario (Fig. 2c andd). For ISBA-CTRIP, negative R 2 coefficients are also reached in north Bolivia, Alaska, north of India and Siberia (south of Lena River). For WGHM, negative R 2 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.

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 (cutoff 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 2 values and reported in Maxwell's colour triangle (Fig. 3 ability (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 pluriannual 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.

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 -1 ), 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 Zam-bezi) Africa, north-east Australia, south Asia (mostly the north of India, especially when using ISBA-CTRIP), northwest 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 landuse 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).

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 climaterelated 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).

Central Amazon Corridor

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. [START_REF] Junk | The Central Amazon floodplain: Ecology of a pulsing system[END_REF][START_REF] Melack | Amazon floodplain hydrology and implications for aquatic conservation[END_REF]. The extension of the flooded area varies from 100 000 to 600 000 km 2 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;[START_REF] Da Silva | Water level dynamics of Amazon wetlands at the watershed scale by satellite altimetry[END_REF], 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;[START_REF] Fan | Global patterns of groundwater table depth[END_REF]. 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 [START_REF] Frappart | The spatio-temporal variability of groundwater storage in the Amazon River Basin[END_REF], which constitutes a large part of the TWS changes observed by GRACE [START_REF] Frappart | The spatio-temporal variability of groundwater storage in the Amazon River Basin[END_REF].

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 andb), 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 -1 ).

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 [START_REF] Frappart | The spatio-temporal variability of groundwater storage in the Amazon River Basin[END_REF]. 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. [START_REF] Decharme | Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications[END_REF][START_REF] Scanlon | Global models underestimate large decadal declining and rising water storage trends relative to GRACE satellite data[END_REF][START_REF] Scanlon | Tracking seasonal fluctuations in land water storage using global models and GRACE satellites[END_REF]. Large floodplains may indeed delay the water transport for several months (e.g. [START_REF] Prigent | Satellite-derived global surface water extent and dynamics over the last 25 years (GIEMS-2)[END_REF], through storage and percolation from the surface towards the aquifer (e.g. [START_REF] Lesack | Flooding hydrology and mixture dynamics of lake water derived from multiple sources in an Amazon floodplain lake[END_REF][START_REF] Bonnet | Floodplain hydrology in an Amazon floodplain lake (Lago Grande de Curuaí)[END_REF][START_REF] Frappart | The spatio-temporal variability of groundwater storage in the Amazon River Basin[END_REF]. Groundwater stores excess water during wet periods and sustains rivers and floodplains during low-water periods (e.g. [START_REF] Lesack | Water balance and hydrologic characteristics of a rain forest catchment in the central Amazon basin[END_REF]. Groundwater systems have also been shown to convey seasonal anomalies (for example, droughts) for several years at local (e.g. [START_REF] Tomasella | The water balance of an Amazonian micro-catchment: The effect of interannual variability of rainfall on hydrological behaviour[END_REF] 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.

Upper São Francisco

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 2 . 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 3 . 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) [START_REF] Lucas | Significant baseflow reduction in the sao francisco river basin[END_REF]. As a result of a prolonged drought lasting from 2002 to 2017 [START_REF] Freitas | Drought assessment in São Francisco river basin, Brazil: characterization through SPI and associated anomalous climate patterns[END_REF], 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 [START_REF] De | Hydroelectric production from Brazil's São Francisco River could cease due to climate change and inter-annual variability[END_REF]. 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 [START_REF] Ndehedehe | Assessing land water storage dynamics over South America[END_REF], allowing the impact of prolonged droughts on the water supply in a vulnerable region to be better quantified (Paredes-Trejo et al., 2021).

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 hy- drological 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.

Zambezi-Okavango

Study area

The Zambezi River basin, located in south tropical Africa, drains an area of 1.4 × 10 6 km 2 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 %) [START_REF] Vörösmarty | Modeling basin-scale hydrology in support of physical climate and global biogeochemical studies: An example using the Zambezi River[END_REF]. 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

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

Congo

Study area

The Congo basin is the second-largest river basin in the world, with a drainage area of ∼ 3.7 × 10 6 km 2 and an average annual discharge of ∼ 40 500 m 3 s -1 [START_REF] Laraque | Recent budget of hydroclimatology and hydrosedimentology of the congo river in central Africa[END_REF]. Despite its importance, the Congo River basin is scarcely studied [START_REF] Alsdorf | Opportunities for hydrologic research in the Congo Basin[END_REF], though a growing interest has arisen over the past decade, substantially due to advances in satellite hydrology (e.g. Papa et al., 2022;[START_REF] Paris | Monitoring Hydrological Variables from Remote Sensing and Modeling in the Congo River Basin[END_REF][START_REF] Schumann | Reviewing Applications of Remote Sensing Techniques to Hydrologic Research in Sub-Saharan Africa, with a Special Focus on the Congo Basin, Congo Basin Hydrology, Climate, and Biogeochemistry: A Foundation for the Future[END_REF]. With an average rainfall around 1500 mm -1 , 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 [START_REF] Alsdorf | Opportunities for hydrologic research in the Congo Basin[END_REF], leading to a bimodal river discharge (Kitambo et al., 2022). The "Cuvette centrale" is a topographic depression lo- cated at the centre of the basin, harbouring wetlands covered by rainforests permanently or periodically flooded [START_REF] Becker | Satellite-based estimates of surface water dynamics in the Congo River Basin[END_REF]. The Congo floodplain hydrodynamics are disconnected from the main river, with much less variability observed throughout the year [START_REF] Alsdorf | Opportunities for hydrologic research in the Congo Basin[END_REF]. 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 [START_REF] Becker | Satellite-based estimates of surface water dynamics in the Congo River Basin[END_REF].

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 [START_REF] Laraque | Recent budget of hydroclimatology and hydrosedimentology of the congo river in central Africa[END_REF], as well as in discharge times series at 7.5 and 13.5 years [START_REF] Labat | Recent advances in wavelet analyses: Part 2-Amazon, Parana, Orinoco and Congo discharges time scale variability[END_REF]. 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 [START_REF] Pfeffer | Analysis of the interannual variability in satellite gravity solutions: detection of climate modes fingerprints in water mass displacements across continents and oceans[END_REF]. 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) [START_REF] Labat | Recent advances in wavelet analyses: Part 2-Amazon, Parana, Orinoco and Congo discharges time scale variability[END_REF]. 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 [START_REF] Laraque | Recent budget of hydroclimatology and hydrosedimentology of the congo river in central Africa[END_REF]; (ii) errors in runoff and evapotranspiration fluxes; or (iii) unresolved underground processes, including preferential flow along faults (Fig. 1 in [START_REF] Garzanti | Congo River sand and the equatorial quartz factory[END_REF].

North Australia

Study area

The climate of north Australia is characterised by a wet season lasting from November to April, subject to intense thun- derstorms and cyclones, with virtually no precipitation during the remainder of the year [START_REF] Smith | Characteristics of the northern Australian rainy season[END_REF]. 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 [START_REF] Petheram | Flow characteristics of rivers in northern Australia: implications for development[END_REF][START_REF] Smerdon | Identifying the contribution of regional groundwater to the baseflow of a tropical river (Daly River, Australia)[END_REF]. Groundwater plays an essential role in north Australia as it sustains rivers and vegetation, through baseflow and water uptake for plant transpiration [START_REF] Lamontagne | Groundwater use by vegetation in a tropical savanna riparian zone (Daly River, Australia)[END_REF][START_REF] O'grady | Groundwater use by riparian vegetation in the wet-dry tropics of northern Australia[END_REF]. Significant interannual variability, principally related to ENSO in the north of the continent, has been observed in rainfall [START_REF] Cai | Fluctuations of the relationship between ENSO and northeast Australian rainfall[END_REF][START_REF] Sharmila | Mechanisms of multiyear variations of Northern Australia wet-season rainfall[END_REF], river discharge [START_REF] Chiew | El Nino/Southern Oscillation and Australian rainfall, streamflow and drought: Links and potential for forecasting[END_REF]Ward et al., 2010) and terrestrial water storage [START_REF] Xie | Multi-climate mode interactions drive hydrological and vegetation responses to hydroclimatic extremes in Australia, Remote Sens[END_REF]. 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 [START_REF] Xie | Multi-climate mode interactions drive hydrological and vegetation responses to hydroclimatic extremes in Australia, Remote Sens[END_REF].

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 [START_REF] Yang | Spatio-temporal evaluation of water storage trends from hydrological models over Australia using GRACE mascon solutions[END_REF], 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 [START_REF] Xie | Multi-climate mode interactions drive hydrological and vegetation responses to hydroclimatic extremes in Australia, Remote Sens[END_REF] because they convey more information about water availability in the soils and aquifers, especially when associated with SMOS (Soil Moisture and Salinity) measurements [START_REF] Li | Global GRACE data assimilation for groundwater and drought monitoring: Advances and challenges[END_REF].

Central United States: Ogallala aquifer

Study area

The Ogallala, or High Plains, aquifer covers a surface area of about 450 000 km 2 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 [START_REF] Houston | Geodatabase Compilation of Hydrogeologic, Remote Sensing, and Water-Budget-Component Data for the High Plains Aquifer[END_REF]. Groundwater abstractions for irrigation began in Texas in the 1930s [START_REF] Luckey | Water-level and saturated thickness changes, predevelopment to 1980[END_REF] and exceeded recharge over much of the central and southern parts of aquifer in the 1950s [START_REF] Luckey | Hydrogeology, water use, and simulation of flow in the High Plains Aquifer in north-western Oklahoma, southeastern Colorado, southwester Kansas, northeastern New Mexico, and northwestern Texas, Water-Resources[END_REF], 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 [START_REF] Haacker | Water level declines in the High Plains Aquifer: Predevelopment to resource senescence[END_REF]. At current depletion rates, a large part of irrigation (about 30 %) may not be supported in the coming decades (Scanlon et al., 2012b;[START_REF] Haacker | Water level declines in the High Plains Aquifer: Predevelopment to resource senescence[END_REF][START_REF] Steward | Peak groundwater depletion in the High Plains Aquifer, projections from 1930 to 2110[END_REF][START_REF] Deines | Transitions from irrigated to dryland agriculture in the Ogallala Aquifer: Land use suitability and regional economic impacts[END_REF].

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 [START_REF] Rateb | Comparison of groundwater storage changes from GRACE satellites with monitoring and modeling of major US aquifers[END_REF], 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 [START_REF] Brookfield | In situ and GRACEbased groundwater observations: Similarities, discrepancies, and evaluation in the High Plains aquifer in Kansas[END_REF] 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 [START_REF] Seyoum | Monitoring and comparison of terrestrial water storage changes in the northern high plains using GRACE and in-situ based integrated hydrologic model estimates[END_REF]. 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.

North of India

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 -1 respectively (e.g. [START_REF] Bhanja | Validation of GRACE based groundwater storage anomaly using in-situ groundwater level measurements in India[END_REF]. The average population density ranges from 26-250 persons per km 2 in the north-west of India to over 1000 persons per km 2 in the north-east of India [START_REF] Dangar | Causes and implications of groundwater depletion in India: A review[END_REF]. India is the largest groundwater user in the world, with an annual withdrawal of 230 km 3 for irrigation, used essentially for rice, wheat, sugarcane, cotton and maize cultures [START_REF] Mishra | Groundwater depletion and associated CO2 emissions in India[END_REF][START_REF] Xie | Multi-climate mode interactions drive hydrological and vegetation responses to hydroclimatic extremes in Australia, Remote Sens[END_REF]. 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 -1 [START_REF] Mishra | Groundwater depletion and associated CO2 emissions in India[END_REF][START_REF] Dangar | Causes and implications of groundwater depletion in India: A review[END_REF]. 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.

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 (-7.71 ± 0.71 mm yr -1 ) 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 [START_REF] Blazquez | Satellite characterization of water mass exchange between ocean and continents at interannual to decadal timescales[END_REF]. Groundwater abstractions were however found to be the dominant driver of water mass losses across northern India (e.g. [START_REF] Xiang | Groundwater storage changes in the Tibetan Plateau and areas revealed from GRACE satellite gravity data[END_REF]. Numerous studies have reported a good agreement between in situ groundwater level measurements and GRACE TWS measurements in the north of India (e.g. [START_REF] Bhanja | Validation of GRACE based groundwater storage anomaly using in-situ groundwater level measurements in India[END_REF][START_REF] Dangar | Causes and implications of groundwater depletion in India: A review[END_REF]. 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 [START_REF] Xie | Multi-climate mode interactions drive hydrological and vegetation responses to hydroclimatic extremes in Australia, Remote Sens[END_REF]. 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. [START_REF] Asoka | Relative contribution of monsoon precipitation and pumping to changes in groundwater storage in India[END_REF][START_REF] Bhanja | Groundwater rejuvenation in parts of India influenced by water-policy change implementation[END_REF].

North of the Black Sea

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 -1 at the north-east of the catchment (Russia) to more than 3000 mm yr -1 at the west (south Austria, Slovenia, Croatia) [START_REF] Rouholahnejad | Water resources of the Black Sea Basin at high spatial and temporal resolution[END_REF](Rouholahnejad et al., , 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 (< -3 • C) in the Krasnodar region (southwest Russia), to over 15 • C at the south of the catchment (north of Turkey) [START_REF] Rouholahnejad | Water resources of the Black Sea Basin at high spatial and temporal resolution[END_REF](Rouholahnejad et al., , 2017)). Land use in the Black Sea catchment is dominated by agriculture [START_REF] Rouholahnejad | Water resources of the Black Sea Basin at high spatial and temporal resolution[END_REF](Rouholahnejad et al., , 2017)).

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 [START_REF] Rouholahnejad | Water resources of the Black Sea Basin at high spatial and temporal resolution[END_REF](Rouholahnejad et al., , 2017)).

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 [START_REF] Werth | Integration of GRACE mass variations into a global hydrological model[END_REF]. 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 [START_REF] Xie | Integrating groundwater irrigation into hydrological simulation of India: Case of improving model representation of anthropogenic water use impact using GRACE[END_REF]. The assimilation of GRACE and GRACE-FO observations into global hydrological models can also increase the model performance with various applications (see [START_REF] Soltani | Review of assimilating GRACE terrestrial water storage data into hydrological models: Advances, challenges and opportunities[END_REF], for a review). Among them, GRACE data assimilation can be used to increase the accuracy (e.g. [START_REF] Zaitchik | Assimilation of GRACE terrestrial water storage data into a land surface model: Results for the Mississippi River basin[END_REF] or resolution (e.g. [START_REF] Kumar | Assimilation of gridded GRACE terrestrial water storage estimates in the North American Land Data Assimilation System[END_REF] 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;[START_REF] Li | Global GRACE data assimilation for groundwater and drought monitoring: Advances and challenges[END_REF]. Significant advances would be expected from the generalisation of such ap-proaches in a dedicated framework (e.g. [START_REF] Condon | Global groundwater modeling and monitoring: Opportunities and challenges[END_REF]Gleeson et al., 2021).

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 1 .

 1 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).

Figure 2 .

 2 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.

Figure 3 .

 3 Figure3. 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.

Figure 4 .

 4 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.

Figure 5 .

 5 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).

Figure 6 .

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

  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.

Figure 7 .

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

Figure 8 .

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

Figure 9 .

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

Figure 10 .

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

Figure 11 .

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

Figure 12 .

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

Figure

  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).

Hydrol. Earth Syst.Sci., 27, 3743-3768, 2023 https://doi.org/10.5194/hess-27-3743-2023

https://doi.org/10.5194/hess-27-3743-2023Hydrol. Earth Syst. Sci., 27, 3743-3768, 2023

Acknowledgements. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation programme (GRACEFUL Synergy grant agreement no. 855677).

Code and data availability. All code and data necessary to validate the research findings have been placed in a public repository at https://doi.org/10.5281/zenodo.7142392 [START_REF] Pfeffer | supplementary-material-for-egusphere-2022-1032 (v0.0)[END_REF].

Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/hess-27-3743-2023-supplement.

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. Author contributions. All authors contributed to the conceptualisation of ideas presented in the manuscript. JP, AlB, BD and SM provided resources necessary to conduct the research findings. JP carried out the formal analysis. AC provided research supervision and funding acquisition. All authors contributed to the investigation of research findings. JP wrote the original draft. All authors contributed to the review and editing of the manuscript.

Competing interests. The contact author has declared that none of the authors has any competing interests.

Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Financial support. This research has been supported by the H2020 European Research Council (grant no. 855677).

Review statement. This paper was edited by Micha Werner and reviewed by three anonymous referees.