Monitoring Groundwater Storage Changes Using the Gravity Recovery and Climate Experiment (GRACE) Satellite Mission: A Review

Abstract

The Gravity Recovery and Climate Experiment (GRACE) satellite mission, which was in operation from March 2002 to June 2017, was the first remote sensing mission to provide temporal variations of Terrestrial Water Storage (TWS), which is the sum of the water masses that were contained in the soil column (i.e., snow, surface water, soil moisture, and groundwater), at a spatial resolution of a few hundred kilometers. As in situ level measurements are generally not sufficiently available for monitoring groundwater changes at the regional-scale, this unique dataset, combined with external information, is widely used to quantify the interannual variations of groundwater storage in the world’s major aquifers. GRACE-based groundwater changes revealed significant aquifer depletion over large regions, such as the Middle East, the northwest India aquifer, the North China Plain aquifer, the Murray-Darling Basin in Australia, the High Plains, and the California Central Valley aquifers in the United States of America (USA), but were also used to estimate groundwater-related parameters such as the specific yield, which relates groundwater level to storage, or to define the indices of groundwater depletion and stress. In this review, the approaches used for estimating groundwater storage variations are presented along with the main applications of GRACE data for groundwater monitoring. Issues that were related to the use of GRACE-based TWS are also addressed.

1. Introduction

Globally, water volume stored in the aquifers is estimated to be about 23.4 × 10⁶ km³, i.e., around half of the fresh water that is only about 2.5% of the total water present on Earth [1]. Aquifer systems play an important role in the global hydrological and the biogeochemical cycles and ecosystems sustainability [2,3,4,5]. Groundwater is a key resource supplying water to billions of people and sustaining agricultural, industrial, and domestic activities [6,7]. It is often the last fresh water resource available for supplying water for domestic use and irrigation after the depletion of surface water in semi-arid areas and densely populated countries [7].

The development of groundwater use has led to depleted levels in many regions around the world [8,9,10]. A global groundwater depletion of 4500 km³ (or 41.4 km³·year⁻¹) was estimated between 1900 and 2008 [11], which represents a major threat to global water security, thus potentially causing a decline in agricultural productivity and energy production. It might also result in conflicts around the world [12]. In the Unites States (US) alone, the groundwater depletion was estimated to be of 1000 km³ (9.2 km³·year⁻¹), with an increase since 1950, and a maximum depletion rate of 24 km³·year⁻¹, observed since 2000 [13]. Globally, the groundwater depletion rate increased from 126 ± 32 km³·year⁻¹ to 283 ± 40 km³·year⁻¹ between 1960 and 2000. In 2000, groundwater depletion reached 39 ± 10% of the global yearly groundwater abstraction [10]. This increase in groundwater use is attributed to growth of the population and an increase of the economic activity, e.g., expansion of the irrigated areas. According to Zekster et al. [14], 20% of the global withdrawals for irrigation, 40% of the total withdrawals for industrial activities, and 50% of the total withdrawals by municipalities come from the groundwater. Total groundwater withdrawals range from 600 to 1100 km³·year⁻¹ (i.e., one-fifth to one third of the total global freshwater abstraction) [14,15,16]. The consumption of groundwater for irrigation is estimated to be 545 km³·year⁻¹ (43% of the total water used for irrigation) [7]. Groundwater depletion (i.e., groundwater abstraction exceeds groundwater recharge), which occurs over a long time-period and extensive areas, has significant impacts on natural streamflow, groundwater level, and, as a consequence, even on land subsidence and salt water intrusion in deltaic areas, in turn, causing the degradation of ecosystems and the deterioration of groundwater quality [10,17]. Globally, groundwater is currently overexploited: its footprint exceeds the actual aquifers area by about 3.5 times. About 1.7 billion of the world’s population lives in areas where groundwater resources are under stress [9].

Climate change is expected to have substantial direct impacts on groundwater water recharge, and hence, on renewable groundwater resources and indirect ones through changes in groundwater use [16,18,19]. Direct impacts relate to climate effects on other reservoirs (droughts and floods intensity and duration in the surface reservoir for instance) and fluxes (rain and snow falls an evapotranspiration) and on land cover. Indirect effects are related to human activity, e.g., expansion of rain-fed and irrigated agriculture [19]. Numerical simulations that were based on climate change scenarios, and taking population density into account, revealed that between 15% and 20% of the world population will be affected by groundwater decreases by the year 2050 [16].

In this context, long-term reliable, compact and frequent observations of the aquifer levels are necessary to monitor groundwater storage and its changes over time. The launch of the Gravity Recovery and Climate Experiment (GRACE) mission in March 2002 offered new perspectives for the monitoring of terrestrial part of the hydrological cycle, including groundwater resources [20]. The GRACE mission measures the tiny variations in the Earth gravity field that are related to changes in Terrestrial Water Storage (TWS) [20,21]. Anomalies of TWS from GRACE are commonly used in hydrological applications from regional to global scales (see [22,23] for recent reviews). Alone or in combination with external datasets, the GRACE data provide a new source of information on extreme climatic events, such as exceptional droughts and floods [24,25,26], on snow accumulation [27], and on basin-scale hydrological fluxes, such as fresh water discharge [28,29] or evapotranspiration [30,31].

Using in situ aquifer level measurements from the High Plains aquifer, Central US, and land surface model outputs, [32] identified the potential of the GRACE mission to detect changes in groundwater (GW) storage using the GRACE data. In their study, they show that the annual amplitude of the GW changes from in situ measurements was higher than the difference between the GRACE-based TWS and the soil moisture (SM) from a land surface model, taking into account the expected accuracy of GRACE and the uncertainty on the model outputs.

In this review, the GRACE mission is described, i.e., on-board instruments, orbital characteristics, and products that were made available to the scientific community. The different approaches that were used to estimate GW changes from GRACE-based TWS are described. Next, the error budget and the validation of the GW storage changes that were estimated using the GRACE data are presented. The main results and main scientific applications are subsequently summarized. The discussion section will present the current limitations of the GRACE-derived GW changes, and which improvements can be expected with the future launches of the GRACE Follow-On (GRACE FO) and the Surface Water and Ocean Topography (SWOT) missions in 2018 and 2021, respectively.

2. The GRACE Mission

2.1. Orbit and On-Board Instruments

The GRACE satellite gravity mission was jointly launched in 17 March 2002 by the American National Aeronautics and Space Administration (NASA) and the German Aerospace Center (Deutsches Zentrum für Luft- und Raumfahrt-DLR). Its scientific mission ended on 27 October 2017. The GRACE mission consists of two co-orbital satellites, which were orbiting at a relatively low average altitude of around 450 km with an inclination of 89.5°. GRACE was unique among all of the the Earth observation satellite missions, since it was composed of two vehicles that followed each other at a distance of about 220 km (Figure 1). The instantaneous measurement of distance and relative speed between the two satellites was performed by a radar telemeter operating at K-band (K-Band microwave Ranging or KBR) with a great accuracy of 1 μm·s⁻¹ [21]. The inter-satellite distance is a function of the differences in gravitational acceleration that was experienced by each GRACE satellite [33]. The payload of each satellite was also composed of a three-axis accelerometer that provides information on the dynamic effects, including non-dissipative or conservative forces, mainly meaning solar and Earth radiation pressure, and atmospheric drag. After removing the non-gravitational effects, variations in distance and inter-satellite distance changing rate were used to estimate variations in the geopotential along the track of the twin satellites. Variations in the Earth’s gravity related to density heterogeneities and topography are reflected in changes in the observed inter-satellite distances along the orbit (Figure 2).

The payload was composed of five instruments on board each satellite component of the GRACE mission [21]:

• the K-band ranging (KBR) system, which provides distance measurements between the two satellites, with an accuracy of 10 μm, using the phases of carrier electromagnetic waves in the K and the Ka bands, at frequencies of 26 and 32 GHz;
• the Ultra-Stable Oscillator (USO), which generates electromagnetic waves in the K-band for the KBR system at the desired frequency;
• the SuperSTAR accelerometers (ACC) accurately measures the non-conservative forces applied to each satellite along three axes;
• the Stellar Camera ASSEMBLY (SCA) determines the orientation of the satellite relatively to the position of fixed stars; and,
• The Black-Jack GPS receivers and Instrument Processing Unit provides the three components of the position and velocity of each of the satellites.

GRACE was the first mission to use the principle of two co-orbital satellites in pursuit, as proposed by [35], allowing for estimating the spatial and temporal variations of the gravitational field (i.e., the geoid), which reflect mass changes in the Earth system over time scales ranging from months to a decade [20]. The GRACE mission that was measured the sum of the effects of all the changes in mass, which were integrated along the vertical axis, of the solid Earth and of redistributions of mass within the superficial fluid envelopes [20].

2.2. Data from the GRACE Mission

Three official processing centers that form the GRACE Science Data Center (SDS) produce the level 1 products and level 2 solutions derived from measurements of the GRACE mission. They are:

• the Center for Space Research (CSR) in Austin, Texas, United States,
• the GeoForschungs Zentrum (GFZ) in Potsdam, Germany, and
• the Jet Propulsion Laboratory (JPL) in Pasadena, California, United States.

Level 1b products consist of already processed positions and velocities that were measured by the on-board GPS, accelerometers, and the accurate K-band measurements of variations in distance between the two vehicles. These measurements are used to compute the monthly gravity field models or level 2 products that were expressed in geoid height or in terms of Equivalent Water Height (EWH, e.g., [36]). These products are made available by the GFZ’s Integrated System Data Center (ISDC—[37]) and JPL’s Physical Oceanography Distributive Active Data Center (PODAAC—[38]).

Other level 2 GRACE solutions can also be obtained by other research groups who might use different approaches at temporal resolutions from one day to one month in the form of spherical harmonic coefficients (global approaches) or spatial grids (local or regional approaches).

The level 2 GRACE products are obtained using a dynamic approach that is based on the Newtonian formulation of the satellite’s equation of motion in an inertial frame, with the Earth’s center as origin, combined with a dedicated modeling of gravitational and non-conservative forces acting on the spacecraft [33]. During the process, the known gravitational effects (changes in atmosphere and oceans masses, pole movement, and tides) are removed from observations using a priori information from meteorological and global ocean circulation models, as well as the non-gravitational forces that were measured by the on-board accelerometers [40,41]. The residual values that form GRACE Stokes coefficients represent, over land, represent the sum of the solid Earth contribution, the continental water storage, the errors from the correcting models, and the noise. More detailed on the phasical processed reflected in GRACE time-variable data can be found in [42]. These Stokes coefficients are first centered to suppress the contribution of the solid Earth, either removing their average that was computed over a period of several or the coefficients from a static geoid model and then converted into equivalent water storage. This processing is described in [36,42].

2.3. Accuracy and Spatial Resolution of GRACE-Based Products over Land

The GRACE mission provided novel information for continental hydrology research: the terrestrial water storage or integrated water content (i.e., the sum of the water contained in the soil column from the following hydrological reservoir: surface water, soil water, groundwater, snow). An early study showed an expected measurement accuracy of a few millimeters of equivalent water height (EWH), in terms of mass surface density (and for a reference water density of 1000 kg.m⁻³), over areas of 400 km by 400 km, which was based on Land Surface Model (LSM) outputs. The TWS retrieval was expected to be affected by the presence of noise in the shorter wavelengths [36]. In addition, spectral truncation errors increase as the area decreases. Based on LSM outputs and the expected accuracy of the GRACE land water solutions, it was estimated that the changes in TWS could be detected if they exceed 1.5 cm of EWH over an area of 200,000 km² [43]. The accuracy of the GRACE land water solutions was expected to be around 0.7 cm EWH for a drainage area of 400,000 km² and 0.3 cm for a drainage area of 4,000,000 km² [44]. Current GRACE products have a spatial resolution of a few hundred kilometers (around 200 km for the Mascons and the regional solutions and 330 km for the releases 03 to 05 for a typical degree of truncation of 60). Errors were estimated to be around 40 mm at the Equator, and decreasing to 15 mm in polar regions [45].

3. Estimating Groundwater Storage Using the GRACE Data: Different Approaches

3.1. The Direct Approach

This straightforward technique consists of removing to the GRACE-based anomaly of TWS (ΔTWS) the contributions of the different hydrological reservoirs, i.e., surface (ΔW_{Surface water}) and soil (ΔW_{Soil water}) waters, snow (ΔW_Snow) and groundwater (ΔW_Groundwater):

ΔTWS = ΔW_{Surface water} + ΔW_{Soil water} + ΔW_Snow + ΔW_Groundwater

(1)

Equation (1) hence becomes:

ΔW_Groundwater = ΔTWS − ΔW_{Surface water} − ΔW_{Soil water} − ΔW_Snow

(2)

The contributions from the individual hydrological reservoirs are obtained from external sources, e.g., hydrological model outputs, in-situ data, or other remotely sensed observations. Figure 3 presents an example of terrestrial trends that were estimated over 2004–2010 from TWS from GRACE with and without the Glacial Isostatic Adjustement (GIA) correction from [46] applied, SM and snow from Global Land Data Assimilation System—National Centers for Environmental Prediction (NCEP), Oregon State University (Department of Atmospheric Sciences), Air Force, Hydrology Lab (GLDAS-NOAH), and the difference between TWS from GRACE corrected from GIA and SM and snow. GIA is the response of the solid Earth to the changes in surface loads due to the melting of large-scale ice sheets and glaciers since the last glacial maximum that occurred ~21,000 years [47]. Over land, the adjustment of the Earth surface to changing loads causes vertical uplifts with the largest trends over Fennoscandia, Northern America, Patagonia, Greenland, and Antarctica. Vertical uplift rates of ~10 mm·year⁻¹ were recorded in Fennoscandia and near the Hudson Bay [47,48].

Several combinations to solve (2) using in situ records, model outputs, and remotely sensed observations are detailed below for a wide range of hydro-climatic conditions.

3.1.1. Regions Where TWS Is Limited to Soil Water Storage

The simplest cases relate to arid and semi-arid environments or in areas where surface waters are negligible and snow inexistent. If a dense soil moisture network exists in areas that are larger than several hundred square kilometers as Illinois or the region of the High Plains aquifer in the United States of America (USA), providing information on the water contained in the soil at several depths, groundwater storage changes can be inferred removing the anomaly of soil water content to the anomaly of TWS that was derived from GRACE [49,50,51]. Most of the times, such networks do not exist. Therefore, hydrological model or land surface model outputs are used instead, such as the ones from GLDAS [52], the WaterGAP Global Hydrology Model (WGHM) [53], and the World-Wide Water Resources Assessment (W3RA) model [54], for instance, for the High Plains aquifer [55], for the Canning aquifer in northwest Australia [56], and within the Middle East [57].

3.1.2. More Complex Environments

In more complex environments, water storage anomalies from the other hydrological reservoirs need to be taken into account using in situ data, hydrological model outputs, or remotely sensed observations, especially when surface water reservoir cannot be neglected. In-situ data were used, in combination with GLDAS outputs for SM, for monitoring the GW changes in the semi-arid Murray Darling Basin, Australia, where all SW is composed of reservoirs, lakes, weirs, and in-channel storage [58], as well as the prone to flood Bengal Basin in Bangladesh [59]. Additional information on snow from models was also taken into account to estimate GW anomalies in California Central Valley, USA [60], and in the Great Lakes Water Basin in Canada and USA [61]. Due to the lack of available in-situ data in remote and complex environments, GW retrievals from GRACE and outputs from hydrological models were performed using; for instance, GLDAS over northern India [62,63] or WGHM in the Niger and the Murray Darling Basins [64,65]. Remote sensing data were also used to provide information on the surface reservoir. Radar altimetry or radar altimetry, coupled with satellite images, was used to remove SW storage variations of lakes to TWS in lakes water basin, such as the East African Lakes [66,67] or the American Great Lakes [61], and in large aquifer regions comprising lakes and reservoirs, such as the Middle East [68,69]. In large drainage basins that were covered with extensive floodplains, surface water storage, which was obtained from the combination of altimetry-based water levels with inundation extent from satellite images, were removed from TWS to retrieve GW changes. This method was applied in the Negro River Basin, which is one of the largest tributary of the Amazon [70], and in the Lower Ob’ Basin in Siberia [71]. In this latter case, the snow storage was also derived from the GRACE data using an inverse approach [72,73].

3.2. Calibration and/or Assimilation into Hydological Models

GRACE data were also incorporated into groundwater, hydrological, and land surface models to improve the simulations of the terrestrial water cycle, including groundwater storage through calibration and data assimilation processes from regional to global scales. Due to the lack of in situ measurements of aquifer levels, hydrological models are imperfect in representing the characteristics of the groundwater changes, which is due to uncertainties of model parameters, limitations in the structure of these models, and uncertainty of input (forcing) data [74]. In spite of the coarse resolution of the GRACE Level-2 products (around 300 km but often re-sampled at 1° or 0.5°), the GRACE-based spatial pattern of GW and associated monthly temporal variations are included in the tuning process of groundwater, hydrological, and land surface models alone or as an additional regularization mechanism. GRACE data were used to constrain the hydraulic head to determine the best values for the model parameters, such as the hydraulic conductivity. Several calibration approaches were applied in various environments—an evolutionary multiobjective optimization on the regional groundwater model that was based on MODFLOW for the Edwards–Trinity Plateau and Pecos Valley aquifers, in Texas, USA [74], two meta-heuristic algorithms (particle swarm optimization and pattern search) on a groundwater flow model for the Ghaen aquifer, Iran [75], Monte-Carlo technique on a large-scale groundwater model that was based on FEFLOW is constructed for the Qaidam Basin, China [76]. These studies demonstrated that the supplementary constrain on storage changes only provided by the GRACE data highly reduced the uncertainty of the model outputs. Calibration and data assimilation of GRACE data were also performed in large-scale hydrological models. A Pareto-based multi-objective calibration approach was applied to WGHM using the GRACE-derived TWS in 28 large drainage basins. If no comparison was performed against in situ data, the calibration resulted in a decrease of the variability of groundwater storage variations at global-scale [77].

The GRACE-based TWS were also assimilated in WGHM through an ensemble-based Kalman filter approach in the Mississippi Basin. Uncertainties in the groundwater model outputs are reduced after one year of simulation, even if they remain large in the region of the High Plains aquifer due to extensive irrigation and groundwater abstraction that are poorly described in the model [78]. Incorporation GRACE-based TWS into land surface models also improved hydrological parameters and fluxes. The calibration of the Community Land Model (CLM) 3.0 [79], coupled with an unconfined aquifer model [80] with a multi-objective approach based on Monte-Carlo simulation framework using TWS from GRACE and base-flow data, improved the simulation of the water table depth. This approach, tested in Illinois, USA, is likely to be extended worldwide [81]. Assimilation of GRACE-based TWS in the Catchment Land Surface Model (CLSM) [82] using an ensemble Kalman smoother method allowed better modeling of groundwater storage variations in the Mississippi River Basin and in 18 river basins in western and central Europe [83,84]. The impact of GRACE error correlation structure on the assimilation of GRACE data was very recently studied by [85,86]. Their results indicate that incorporating GRACE data’s full covariance matrices results in more realistic groundwater estimations. Recently, Khaki et al. [86,87] applied GRACE data, through various extensions of assimilation filters and applying covariance localization, and [88] used GRACE and soil moisture data simultaneously in an ensemble-based assimilation framework to update storage estimation of a hydrological model in Australia, while [57,89] assimilate GRACE data into W3RA to estimate the groundwater changes in Bangladesh and the Middle East, respectively.

4. Error Budget and Validation of the GRACE-Based Groundwater Storage

4.1. Error Budget of the GRACE-Based Groundwater Storage

Using Equation (2), the error on GW storage variations (σ_Groundwater) can be written as:

$${\mathsf{\sigma }}_{\mathrm{Groundwater}}=\sqrt{{\left({\mathsf{\sigma }}_{\mathrm{TWS}}\right)}^2+{\left({\mathsf{\sigma }}_{\mathrm{Surface}\mathrm{water}}\right)}^2+{\left({\mathsf{\sigma }}_{\mathrm{Soil}\mathrm{water}}\right)}^2+{\left({\mathsf{\sigma }}_{\mathrm{Snow}}\right)}^2}$$

(3)

where σ_TWS, σ_{Surface water}, σ_{Soil water}, and σ_Snow represent the error on the TWS, surface water, soil water, and snow components, respectively.

Or similarly, on the trend:

$${\mathsf{\sigma }}_{\mathrm{aGroundwater}}=\sqrt{{\left({\mathsf{\sigma }}_{\mathrm{aTWS}}\right)}^2+{\left({\mathsf{\sigma }}_{\mathrm{aSurface}\mathrm{water}}\right)}^2+{\left({\mathsf{\sigma }}_{\mathrm{aSoil}\mathrm{water}}\right)}^2+{\left({\mathsf{\sigma }}_{\mathrm{aSnow}}\right)}^2}$$

(4)

where σ_aTWS, σ_{aSurface water}, σ_{aSoil water}, and σ_aSnow represent the one-sigma trend error on the TWS, surface water, soil water, and snow components, respectively.

Errors of GRACE data are originated from measurements, spatial and spectral leakages, post-processing, and glacial isostatic adjustment (GIA) in areas where this phenomenon is significant. Errors on trend estimates of the GRACE-based GW storage were found to be equal to 5.3 mm·year⁻¹ for a trend of −55.9 mm·year⁻¹ from April 2006 to March 2010 in the California Central Valley aquifer [92], of 0.3 cm·year⁻¹ for a trend of −2.2 cm·year⁻¹ from 2003 to 2010 in the north of China [93], 0.4  km³·year⁻¹ for a trend of −5.6  km³·year⁻¹ from December 2004 to November 2013 in the Colorado Basin [94], or 0.6 mm·year⁻¹ for a trend of −27.2 mm·year⁻¹ from 2003 to 2009 in the Tigris and the Euphrates Basin [68], when computing the error budget presented in Equation (4).

4.2. Validation of the GRACE-Based Groundwater Storage Variations

Validation of the GRACE-based GW storage was performed against in situ well measurements when available. Two different approaches were used: the validation of GRACE-based TWS comparing them to the sum of the soil moisture storage (from in-situ data or from model outputs) to the groundwater storage (from wells) or the direct validation of the GRACE-based groundwater storage anomalies that were obtained using Equation (2) against the well records. These approaches required two preliminary steps. The well depths (h) are proportional to groundwater storage through the specific yield (S_y):

W_Groundwater = S_y h

(5)

Once the conversion is performed, in situ SM and GW storages are either interpolated or are averaged to be compared with GRACE-based TWS or GW storage.

The first approach was applied in Illinois using 16 wells and 16 soil moisture stations from August 2002 to November 2005 on area of 200,000 km² [49], and 16 wells and 19 soil moisture stations from February 2003 to December 2005 on area of 280,000 km² [95] using the GRACE Release-03 data from CSR. Very similar results were found (R = 0.83 in the first study and R = 0.95 and RMSE = 20.3 mm in the second study). The differences that were observed between these two studies can be attributed to difference in observation period, validation datasets, extent of the study area, and data processing (SM were considered to 2 and 1 m depths, respectively). The better results that were obtained are more likely to be accounted by the choice of the period (first GRACE data from 2002 are of lower quality than the ones obtained after February 2003) and a larger study area. Similar validations were performed in the High Plains aquifer (450,000 km² of area) using (i) GRACE Release 03 data from CSR, SM outputs from NOAH LSM [96] and 2700 wells records between January 2003 and December 2005 [55]; (ii) GRACE Release 04 data from CSR, SM outputs from NOAH LSM, and from 78 in-situ stations and 1000 wells records between January 2003 and December 2006 [50]. Correlation of 0.82 and RMSE of 33 mm were found in the first study and R of 0.95 and 0.96 and RMSE of 36 and 38 mm using simulated and in situ SM, respectively, in the second study. The increase of correlation is more likely resulting from the use of a more recent GRACE dataset.

The second approach was applied by the same groups in Illinois and in the High Plains aquifer. The following results were obtained: R = 0.63 in Illinois [49], R = 0.58 [50], R = 0.72, and R = 0.73 using simulated and in situ SM [55] for the High Plains aquifer. In all these cases, R values are lower for GW storage than for TWS. It was also used in several other areas where a large number of in situ data was available. R² equals 0.89 was found between GRACE-based TWS from the average of CSR, GFZ and JPL Release 04 solutions minus SM from GLDAS/NOAH and GW storage that were derived from 935 wells records between January and December 2008 in the state of Gujarat (196,030 km² of area), India [97], and it equals 0.65 and 0.60 (with RMSE of 85 and 49 mm) on storage variations and 0.91 and 0.75 (with RMSE of 23 and 14 mm) on yearly smoothed storage variations between GRACE-based TWS from CSR Release 05 solutions minus SM from the average of the four LSM used in GLDAS and Climate Prediction Center (CPC) LSM [98] and GW storage derived from 87 wells records between January 2005 and December 2010 in Piedmont Plain (54,000 km² of area) and East Central Plain (86,000 km² of area), China [99]. Very good agreement was also found in the Bengal Basin when comparing GW storage obtained removing SW from 298 water level gauge stations, SM from the average of three LSM used in GLDAS from TWS from GRACE CSR Release 04 and GRGS Release 02 solutions and GW from 236 in situ borehole records from January 2003 to December 2007 (R = 0.93 and 0.77, respectively) [59]. Comparisons at annual time-scale were also performed between GW storage that was obtained as the difference between GW storage estimated as the difference between GRACE-based TWS from GRGS Release 01 and from CSR, GFZ, and JPL RL04, surface water storage from in situ measurements, SM from NOAH and in situ measurements from 1462 bores between 2002 and 2007 over the Murray Darling Basin. Correlation of 0.92 was found using the GRGS product [58], and the maximum deviation lower than 5 km³ for variations around 60 km³ using CSR, GFZ, and JPL products [100].

When few in situ aquifer levels were available as in the Amazon Basin, comparisons were performed with surface water levels in areas where the groundwater table permanently reaches the surface. It is the case in the Negro River Basin, one of the largest tributaries of the Amazon, with a drainage area of 700,000 km². In this basin, GW storage that was obtained by the difference between GRACE-based TWS average of CSR, GFZ, and JPL, SW storage that was derived from the multi-satellite observations [101], and SM from the average of the outputs from Land Dynamics (LaD) [102] and WGHM between January 2003 and December 2004, were compared to in situ wells levels and surface water levels derived from radar altimetry. Correlations of 0.76 and 0.73 were found for two regions where the groundwater table permanently reaches the surface [70] (Figure 4).

5. Results

5.1. Groundwater Depletion

The new information on TWS anomaly provided, for the very first time at global scale, by the GRACE mission when combined with external information (see Section 3), allowed for the detection of changes in GWS in the major aquifer worldwide. In most of the major aquifers that were located in arid and semi-arid areas at mid-latitude, large depletion rates were identified, resulting from an excess of water extraction due to intensive pumping when compared with the natural recharge from infiltration. Examples of GWS time-series are presented in Figure 5. Large aquifer depletion rates, down to ~−60 mm·year⁻¹ of GWS loss in EWH, were estimated in North Africa, Middle East, India, China, Australia, and the USA (

Table 1. Annual rates of groundwater depletion in the major aquifers located in arid and semi-arid regions based on the decomposition of GRACE-based TWS.

Region	Area 106 km2	Countries	Estimated Depletion Rate		GRACE Product and Release	SM Product	SW Product	Snow Product	Period (MM/YYYY–MM/YYYY)	Ref
			mm·year−1	km3·year−1
Northwestern Sahara Aquifer System	1	Algeria, Lybia, Tunisia	−0.54  ±  1.40	−0.54  ±  1.40	NA	GLDAS	None	None	01/2003–12/2007	[105]
			−0.81 ± 0.16	−0.81 ± 0.16	CSR 05	None	None	None	01/2003–09/2012	[106]
			−4.48	−4.48	Regional	None	None	None	06/2003–12/2012	[67]
			2.69  ±  0.8	−2.69  ±  0.8	CSR 05	GLDAS	CLM 0.4	None	01/2003–12/2013	[107]
Nubian Sandstone Aquifer System	2.2	Chad, Egypt, Lybia, Sudan	−3.72 ± 0.27	−8.18 ± 0.59	CSR 05	None	None	None	01/2003–09/2012	[106]
			−2.76 ± 0.86	−6.08  ±  1.9	CSR 05	GLDAS	CLM 0.4	None	01/2003–12/2013	[107]
Tigris and Euphrates River Basin	0.75	Turkey, Syria, Iraq, Iran	−27.2 ± 0.6	−20.4 ± 0.45	CSR 05	GLDAS	Altimetry	None	01/2003–12/2009	[68]
Middle East	NA	Irak, Iran, Saudi Arabia, Turkey	NA	−43 ± 3	CSR 05	CLM 4.5	CLM 4.5 + Altimetry	None	02/2003–12/2012	[108]
Northwest India	0.56	India	−40 ± 10	−17.7 ± 4.5	CSR 04	GLDAS	None	GLDAS	08/2002–10/2008	[63]
			−47.7 ± 12	−20.4 ± 7.1	CSR 05	GLDAS	None	GLDAS	01/2003–12/2012	[109]
Norhern India	2.7	India	−20 ± 3	−54 ± 9	CSR 04	CLM	CLM	GLDAS	04/2002–06/2008	[62]
North China	0.37	China	−22 ± 3	−8.3 ± 1.1	CSR 05	GLDAS	None	None	01/2003–12/2010	[93]
Hai River		China	−17.0 ± 4.3	−5.5 ± 1.4	GRGS 2	In situ	In situ	None	01/2003–12/2012	[110]
			−8.3 ± 4.5	−2.7 ± 5.5	CSR 05
Piedmont Plain	0.054	China	−46.5 ± 6.8	−2.5 ± 0.4	CSR 05	GLDAS	None	None	01/2003–07/2013	[99]
East Central Plain	0.086	China	−16.9 ± 1.9	−1.5 ± 0.2	CSR 05	GLDAS	None	None	01/2003–07/2013	[99]
Canning	0.43	Australia	−25.6	−11	GRGS 2	NOAH	None	None	01/2003–12/2009	[56]
Murray Darling Basin	1	Australia	−17 ± 7	17 ± 7	GRGS 1	NOAH	In situ	None	01/2003–12/2007	[58]
			−17.2 ± 4.7	17.2 ± 4.7	CSR 05	WGHM	WGHM	None	01/2003–12/2012	[65]
High Plains Aquifer	0.45	USA	−25.1 ± 2.1	−11.4 ± 0.4	CSR 05	NOAH	In situ	None	03/2003–02/2013	[111]
Sacramento and San Joaquin River Basins	0.154	USA	−20.4 ± 3.9	−3.1 ± 0.6	CSR 04	GLDAS	In situ	NORHSC	10/2002–03/2010	[60]
			−55.9 ± 5.3	−8.6 ± 0.8	GRGS 2	GLDAS	In situ	SNODAS	04/2006–03/2010	[92]
			−44.9 ± 8.5	−6.9 ± 1.3	CSR 04	GLDAS	In situ	SNODAS	04/2006–03/2010	[92]
Colorado River Basin	0.64	USA	8.75 ± 0.63	−5.6 ± 0.4	CSR 05	GLDAS	In situ	SNODAS	12/2004–11/2013	[94]

). In many of these regions, high agricultural productivity is obtained through the use of GW for irrigation [12].

As it can be seen in Table 1, large variations can be observed over the same area from different studies. These apparent discrepancies can be accounted for several factors that include the accuracy of the GRACE datasets (processing center, release, global, or regional solutions), the completeness of the external information, the period on which the trends were estimated. The two former reasons will be analyzed in Section 6. Discussion. The latter one is perfectly illustrated in the case of the Sacramento and San Joaquin River Basins that encompass the Californian Central Valley aquifer. Using almost the same datasets, except for the snow, changes in GWS of −3.1 ± 0.6 and −6.9 ± 1.3 km³·year⁻¹ were found over the periods October 2002–March 2010 and April 2006–March 2003, respectively [60,92]. They result from an acceleration of the depletion that started in the mid-2000’s due to the over abstraction that was caused by intensive irrigation practices [103] Similar accelerations in the decline of GWS were also observed in several locations that were related to persistent drought as in Iran [104].

5.2. Determining Groundwater Related Parameters Using GRACE-Based TWS

GRACE-derived TWS also provides useful information on groundwater related parameters. GW storage estimates obtained subtracting SM from the North America Land Data Assimilation System (NLDAS) to GRACE-based TWS, and well measurements were combined in a robust optimization approach, minimizing the negative impact of uncertainty on the different datasets through the incorporation of boundary conditions, that allows for the retrieval of S_y. This method was applied in the interconnected Edward Trinity Plateau and Pecos Valley aquifers (area of 115,000 km²), in Texas, USA, between August 2002 and December, providing consistent S_y values in the different hydrogeological regions [112]. GRACE-based TWS were also used for predicting GW level change based on the downscaling the GRACE data using artificial neural network models. Test performed for different wells that are located across the US demonstrate the potential of this approach to predict water levels at monthly and season time-scales even in the case of extreme climate events [113].

GRACE-based TWS were also used to characterize the impact of both climate variability and man-induced effects on GWS. GRACE-based GWS anomalies were used to define a GRACE GW Index (GGWI) to quantify the aquifer depletion. This index was defined as the normalized deviation in GWS (i.e., the GWS deviation, given by the difference between the monthly value and its climatology, minus its mean and divided by its standard deviation). Applied to the Central Valley aquifer, California, it shows a good capacity to properly identified GW-related droughts when classical indexes, such as Palmer Drought Index (PDSI) or Standardized Precipitation Index (SPI), detect them with an advance of several months [114]. As GW is a threatened resource, GRACE-based TWS were used to assess groundwater through the Renewable Groundwater Stress (RGS) defined as the ratio between groundwater use and groundwater availability. These two quantities were obtained from the trend in GWS anomalies from GRACE (computed as in [94]) and the natural groundwater recharge (estimated as the vertical flux between the aquifer and bottom soil layer of CLM 4.0). This approach was applied to the world’s 37 largest aquifers to classify them into stress regimes are defined as Overstressed, Variable Stress, Human-dominated Stress, and Unstressed depending on the RGS ratio [115]. Similarly, the Total Groundwater Stress (TGS) was defined as the ratio between the volume of water of the aquifer and the trend in GWS anomalies from GRACE. The former term is estimated as the product between the effective thickness of the aquifer, the porosity, and the surface of the aquifer. TGS was used to determine the number of years until the water volume stored in the aquifers is depleted from 25% to 90% of its total capacity [107].

6. Discussion

The GRACE mission provided estimates of TWS anomalies from March 2002 to June 2017 that were offering an unprecedented source of information to study the terrestrial branch of the hydrological cycle. This unique dataset was used to estimate GWS using the approaches presented in Section 3. However, these methods are suffering from limitations due to the errors on the GRACE data, the non-consideration of some hydrological compartments, fluxes, and abstractions.

6.1. Sources of Errors in the GRACE Solutions

These are predominently of three types:

• Aliasing of short periods from correcting a priori models

Short time periods from hours to days of atmosphere and ocean tides are removed using dealiasing techniques based on corrections from models. While atmosphere pressure fields from ECMWF allow for a reasonable de-aliasing of high-frequency changes that are caused by non tidal atmospheric mass redistributions, errors due to tide model are present in the GRACE solutions, especially for diurnal (S1) and semi-diurnal (S2) tides [116,117,118,119]. Aliasing of the S2 tide has a strong impact on the determination on the degree-two order-zero harmonic coefficient of the geo-potential [120].

• North-South striping

The global Level-2 solutions describe the gravity signatures of un-modeled surface mass variations, mainly the seasonal water mass redistributions on ocean and land areas. Unfortunately, they are affected by north-south striping (see Figure 6), which is particularly visible in the tropical band where the coverage of the satellite tracks is poor because of (1) sparse GRACE track sampling in the longitudinal direction due to the polar orbit of the satellites, (2) propagation of systematic errors from the correcting model accelerations [116,117,118], and (3) strong numerical harmonic coefficient correlations that are generated by solving the undetermined systems of normal equations for high-degree Stokes coefficients [90].

• Spectral representation

Representation in spherical harmonics favors frequency description instead of spatial localization [122]. Stokes coefficients that were determined from GRACE measurements are provided as monthly (or 10-days) lists of harmonic coefficients of the geopotential, which were developed up to a fixed degree, such as 60 (equivalent spatial resolution of 333 km). Gibb’s effects create oscillations by omitting larger degree coefficients in the gravity signal reconstruction. This spectral truncation (or loss of energy) does not give access to the details of water mass fluxes at sub-continental scales (i.e., “leakage out”). Moreover, strong variations of hydrological signals in large drainage basins propagate on the entire sphere and pollute regional estimates in surrounding areas by unrealistic oscillations (i.e., “leakage in”) (e.g., see [100,123] for the cases of several large river basins).

To minimize the effects of the two latter sources of error, strategies to restore the GRACE signal in global solutions and local approaches were developed to be an alternative to the issues that are related to the classical global representation in spherical harmonics through the choice a suitable geometry of surface tiles. For the global solutions, the strategies can be directly based on water storage variations from hydrological model or LSM outputs (i.e., additive [124,125] and scaling factor [45] approaches) or indirectly (i.e., multiplicative approach [126], constrained [127] and unconstrained [72] forward modeling). Using a constrained forward modeling approach, the long term changes in GWS were reassessed in the northwest India aquifer [127]. A depletion rate of 31 ± 1 mm·year⁻¹ (14.0 ± 0.4 km³·year⁻¹) was found between January 2005 and December 2010, which is consistent with in situ measurements (28 mm·year⁻¹ or 12.3 km³·year⁻¹), but lower than the previously published ones (see Table 1), revealing a possible overestimate using classical techniques.

The local approaches are kown as Mascons (for mass concentrations) and regional approaches. In the Mascons approach, the variations of water mass are determined in rectangular elements, which are based on radar telemetry KBR, along the satellite tracks of GRACE, above the study area. The adjustment method of “regional” Stokes coefficients uses spatial and temporal constraints that induce a smoothing effect. Monthly Mascons solutions at 0.5° of spatial resolution were made available by CSR [128] and JPL [129]. Reduction of the leakage effect by 50% was obtained. The regional approach adjusts the surface densities of mass based on Level-1 GRACE data. This method is based on precise measurements on the K-band and the variations of speed between the two satellites. Contrary to the Mascons, the representation mode in spherical harmonics is not used. Therefore, the regional solutions are not affected by spectral truncation effects, and, besides, they provide a better geographical localization of hydrological patterns [130]. These solutions are not yet available globally but over continents, such as Africa [67], Australia [131], and South America [132]. Examples of TWS using regional and global approaches are presented in Figure 7.

6.2. Impacts of the Non Inclusion of Hydrological Water Compartments and Fluxes

The direct approach that is commonly used for estimating GWS anomalies assumes to take into account the water volume that is stored in the individual hydrological compartments. In many studies, LSM outputs were used to represent the other hydrological reservoirs than the GW. If the water that is contained in the other different hydrological reservoirs is occasionally accounted for, if available, using in situ data, abstractions for irrigation or the surface water stored in wetlands and on floodplains are generally not taken into account. In agricultural regions, taking into account the effect of irrigation on GWS changes from direct measurements of irrigation (IRR and SM+IRR) or through proportionality of the difference between potential evapotranspiration (PET) and evapotranspiration (ET) rates (Soil Moisture Deficit–SMD) strongly improves the quality of the estimates when compared to the classical approach. In three different regions of the High Plains Aquifer, r² between in situ and GRACE-based estimates increased from 0.18 to 0.73 and 0.79 in the Northern High Plains Aquifer, 0 to 0.81 and 0.93 in the Central High Plains Aquifer and from 0.02 to 0.44 and 0.72 in the Southern High Plains Aquifer, USA, between 2003 and 2013, taking into account SMS, SMS+IRR, or SMD (both methods give the same results) and IRR [111]. In large river basins with extensive floodplains, SWS cannot be neglected, representing, for instance, between 40% and 50% of the annual variations of TWS as in the Amazon [133,134], Orinoco [135], and Ganges [136,137] basins. Regional SWS water changes can be estimated using multi-satellite observations, such as the synergy between inundation extent from satellite images and either satellite altimetry or Global Digital Elevation Models (GDEM). The first approach was used over the Negro River [70] and the Lower Ob’ [71] basins. It will be soon extended with the future launch of the NASA/CNES Surface Water and Ocean Topography (SWOT) that will provide global maps of water levels at 100 m of spatial resolution over land (in the water mask) [138].

7. Conclusions

During its sixteen years of operation, from March 2002 to June 2017, the GRACE mission provided a novel source of information on TWS variations at unprecedented spatial and temporal resolutions. It offered an exceptional dataset for studying large-scale water mass redistributions, and for the very first time, the opportunity to monitor groundwater changes from regional to global scales. GRACE data were widely used for characterizing the aquifers depletion in vast regions, large drainage areas, and transboundary aquifers, retrieving ground-water related parameters, such as the specific yield and developing new indices to identify aquifer depletion and GW stress. When the GRACE mission ended in 2017, a partnership between NASA and the German Research Centre for Geosciences (GFZ) decided to schedule the launch of the GRACE Follow-On (GRACE-FO) in 2018 in order to ensure continuity of GRACE-type space gravimetry. By using the same twin satellite configuration of low-Earth and the nearly polar orbit of 300–500 km, the GRACE-FO mission will be following its successful predecessor. Additionally, it will carry a demonstrator of laser system to measure the inter-satellite distance and velocity and for an improved precision. It is expected to provide new perspectives in hydrology studies, e.g., refined long-term mass balance estimates of surface water storage and ice sheets. Regional/local approaches have already demonstrated the possible access to spatial (better localization of structures by construction) and temporal (through daily updates using Kalman filter strategies [139,140]) scales that were lower than those that were offered by global solutions. These new data that will offer continuity with the previous GRACE measurements will be very useful for the purpose of global hydrology of global hydrological model calibration and to constrain their operations through assimilation techniques. In spite of its value, GRACE-FO data will still suffer the same limitations, as pointed out by [141]: (i) one-dimensional information for a three-dimensional aquifer; (ii) a coarse resolution incompatible with GW resources management, and (iii) inability to address issues that are related to GW pumping. The synergy with Interferometry Synthetic Aperture Radar (InSAR) data, which provides land subsidence, however, may give access to high resolution and quantitative GW changes maps that could be used to estimate aquifer reaction to pumping, to discriminate between climatic and anthopogenic storages changes and to assess aquifer confinement and dynamics [142]. A first demonstration of the synergy between the two techniques was performed in Central Mexico. Using ALOS-PALSAR images that were acquired between 2007 and 2011, areas subject to ground motion related to GW abstractions were detected. Due to the coarse resolution data, GW losses were only partly observed [142]. This encouraging first result offers new directions for GW studies.