RECOG RL01: correcting GRACE total water storage estimates for global lakes/reservoirs and earthquakes

Observations of changes in terrestrial water storage (TWS) obtained from the satellite mission GRACE (Gravity Recovery and Climate Experiment) have frequently been used for water cycle studies and for the improvement of hydrological models by means of calibration and data assimilation. However, due to a low spatial resolution of the gravity field models, spatially localized water storage changes, such as those occurring in lakes and reservoirs, cannot properly be represented in the GRACE estimates. As surface storage changes can represent a large part of total water storage, this leads to leakage effects and results in surface water signals becoming erroneously assimilated into other water storage compartments of neighbouring model grid cells. As a consequence, a simple mass balance at grid/regional scale is not sufficient to deconvolve the impact of surface water on TWS. Furthermore, non-hydrology-related phenomena contained in the GRACE time series, such as the mass redistribution caused by major earthquakes, hamper the use of GRACE for hydrological studies in affected regions.

In this paper, we present the first release (RL01) of the global correction product RECOG (REgional COrrections for GRACE), which accounts for both the surface water (lakes and reservoirs, RECOG-LR) and earthquake effects (RECOG-EQ). RECOG-LR is computed from forward-modelling surface water volume estimates derived from satellite altimetry and (optical) remote sensing and allows both a removal of these signals from GRACE and a relocation of the mass change to its origin within the outline of the lakes/reservoirs. The earthquake correction, RECOG-EQ, includes both the co-seismic and post-seismic signals of two major earthquakes with magnitudes above M_w9.

We discuss that applying the correction dataset (1) reduces the GRACE signal variability by up to 75 % around major lakes and explains a large part of GRACE seasonal variations and trends, (2) avoids the introduction of spurious trends caused by leakage signals of nearby lakes when calibrating/assimilating hydrological models with GRACE, and (3) enables a clearer detection of hydrological droughts in areas affected by earthquakes. A first validation of the corrected GRACE time series using GPS-derived vertical station displacements shows a consistent improvement of the fit between GRACE and GNSS after applying the correction. Data are made available on an open-access basis via the Pangaea database (RECOG-LR: Deggim et al., 2020a, https://doi.org/10.1594/PANGAEA.921851; RECOG-EQ: Gerdener et al., 2020b, https://doi.org/10.1594/PANGAEA.921923).

1 Introduction

The dynamic global water cycle influences our everyday lives by affecting freshwater availability, weather/climate fluctuations and trends, seasonal variations, anthropogenic water use, and single extreme events such as floods and droughts. Understanding how water is transiently stored and exchanged among the different compartments (groundwater, surface water, soil moisture, etc.) with the help of hydrological models is, therefore, of major societal importance. However, large model uncertainties caused by errors in climate forcings and an incomplete realism of process representations limit the models' ability to accurately simulate water storages and fluxes, making independent observations indispensable for model validation/calibration and data assimilation.

Since 2002 measurements of time-variable gravity obtained from the twin-satellite mission GRACE (Gravity Recovery and Climate Experiment; Tapley et al., 2004) and its successor mission GRACE-Follow-On (GRACE-FO; Flechtner et al., 2016; Kornfeld et al., 2019) have allowed for the determination of column-integrated terrestrial water storage (TWS) changes on a global scale with uniform data coverage (e.g. Pail et al., 2015). However, several challenges are involved with using GRACE for improving hydrological models, among them (1) the low spatial resolution of GRACE, integrating spatially over regions as large as ∼200 000 km² and hampering the representation of concentrated and sub-scale water storage changes, and (2) the fact that gravity observations also contain non-hydrology-related mass variations.

The first problem is caused by the GRACE orbit configuration in combination with unmodelled short-periodic mass changes, resulting in the gravity field models being strongly corrupted by spatially correlated noise. The necessary spatial filtering approach (e.g. Kusche, 2007) inevitably leads to signal loss and to leakage effects resulting in a rather coarse spatial resolution of the gravity field models of a few hundred kilometres.

This limits the investigation of mass variations to rather large-scale processes (Longuevergne et al., 2013), even though small-scale mass variations, whose typical size is smaller than GRACE resolution but large enough in magnitude, can have a strong influence on the total mass change signal (Frappart et al., 2012). Examples are human-controlled reservoirs or natural lakes with strong (seasonal) variations and/or trends. Even though GRACE can “see” these mass changes, they do not necessarily appear exactly at the location of their origin and with the correct magnitude. Thus they can distort the water storage estimate for neighbouring areas or the average over a river basin, as shown in Fig. 1.

Figure 1 Overview of the lake leakage problem with localized changes in water level of the lake/reservoir influencing the estimated water storage in surrounding areas.

This sub-scale mass variability impacts GRACE amplitudes up to 20 % averaged over basins as large as ∼ 200 000 km² (Longuevergne et al., 2013; Farinotti et al., 2015). Although the issue of concentrated hydrological mass variations is not limited to surface waterbodies (e.g. Castellazzi et al., 2018), dam operations and impoundment have a large impact on the water cycle and on continent–ocean exchanges (Chao et al., 2008). Therefore, several publications have focused on removing the impact of surface waterbodies from GRACE total water storage changes. For example, Grippa et al. (2011) removed the influence of surface water storage in the Niger River (derived from altimetry and remotely sensed surface water extent) from GRACE TWS estimates to better be able to compare them to hydrological model output. Tseng et al. (2016) determined mass changes in two Tibetan lakes, combining altimetry, remote sensing, and GRACE estimates, and Zhang et al. (2017) estimated water volume changes for 96 % of the lake area on the Tibetan Plateau by combining an average of ICESat-derived elevation changes from a number of larger lakes with Landsat lake area changes for smaller lakes. Ni et al. (2017) removed the leakage error in GRACE estimates over Lake Volta in Ghana using constrained forward modelling. All these studies represent regional test cases, but a global assessment of the influence of surface waterbody mass change on GRACE data is missing.

Using GRACE data for the evaluation of (global) hydrological models or for combining models and observations by model calibration (Werth and Güntner 2010) and data assimilation (C/DA; Zaitchik et al., 2008; Eicker et al., 2014) without accounting for localized surface water storage can lead to two different kinds of errors. (i) Many global hydrological models do not include a surface water storage compartment at all (Scanlon et al., 2017), and assimilating GRACE TWS into a model that does not explicitly include surface waters will inevitably result in other storage compartments (such as soil moisture or groundwater) becoming distorted by absorbing the observed surface water mass change. Even if a model such as WaterGAP (Müller Schmied et al., 2014) does include a surface water compartment, it might not represent the realistic behaviour of, for example, human-operated reservoirs, and it might be preferable to exclude the reservoir storage from the assimilation. (ii) The leakage effect of localized surface waterbodies might cause an assimilation of the surface water mass change into neighbouring grid cells that should not be affected by it. To our knowledge, no investigation so far has studied the effects of surface waterbodies on GRACE model calibration or data assimilation and how they can best be handled in order to not distort the C/DA results. Having a global correction dataset to clear GRACE water mass changes of the influence of large surface waterbodies (here, lakes and reservoirs) will be immensely helpful for making GRACE estimates more consistent with model output.

Today, extensive information on surface water variations is available from satellite remote sensing. For almost 30 years, satellite altimetry has been providing water levels of large and medium lakes and reservoirs on a global scale (e.g. Birkett, 1995; Berry et al., 2005; Göttl el al., 2016). Several databases make these time series freely available for hydrological applications, among them the Database for Hydrological Time Series of Inland Waters (DAHITI; Schwatke et al., 2015). In addition, optical satellite images are used to derive surface extent of lakes and reservoirs (e.g. Pekel et al., 2016; Klein et al.,2017; Schwatke et al., 2019). Time series from optical sensors can reach a length of up to almost 40 years with a spatial resolution of 250 (MODIS, since 1999), 30 (Landsat, since 1982), and 10 m (Sentinel-2, since 2015) as well as a high temporal resolution with revisit time from 14 (Landsat), over 5 (Sentinel-2), and up to 1 d (MODIS). By combining height and surface area information, time series of storage changes can be derived purely based on remote sensing data (e.g. Schwatke et al., 2020; Busker et al., 2019; Crétaux et al., 2011).

To account for the second challenge in using GRACE data for hydrological studies, namely the removal of all non-hydrology-related mass variations, some effects are typically subtracted using geophysical models either during the computation of the gravity field solutions (e.g. Earth tides, ocean tides, and oceanic/atmospheric mass variations) or in post-processing (e.g. glacial isostatic adjustment (GIA)). However, in addition to this, the mass redistribution caused by the crustal deformation following large earthquakes is also contained in the GRACE observations masking hydrological phenomena in the affected regions. Several studies have highlighted GRACE's usefulness for estimating large earthquakes (M_w>9.0; e.g. Panet et al., 2007; Broerse, 2014; Einarsson et al., 2010; Einarsson, 2011; L. Wang et al., 2012) and have also identified co- and post-seismic earthquake signals in GRACE data with a lower magnitude (e.g. Han et al., 2016; Zhang et al., 2016), down to magnitude 8.3 (Chao and Liau, 2019). For example, the onset of the Sumatra–Andaman earthquake in December 2004 (magnitude 9.1) was analysed using, among others, differences of monthly gravity solutions (Han et al., 2006), wavelet analysis (Panet et al., 2007), Bayesian approaches (Einarsson et al., 2010), or normal modes (Cambiotti et al., 2011). At time of writing, the German GeoForschungsZentrum in Potsdam (GFZ) is the only processing centre that provides a total water storage (level 3) dataset corrected for earthquakes (Boergens et al., 2019); however, a data-based global earthquake correction for different GRACE solutions is not available yet.

To account for both the localized surface water storage in lakes/reservoirs and the earthquake signal, we present the first release of a new global correction dataset, RECOG (REgional COrrections for GRACE) RL01, which can be used for disaggregation of the integral GRACE water storage estimates in addition to applying standard corrections such as GIA models and the atmosphere/ocean de-aliasing products. The term RECOG refers to the fact that all effects included in the data product are localized phenomena that nevertheless influence a larger region around them. The surface water correction (RECOG-LR) was computed from forward-modelling altimetry and remote sensing observations and can be used (a) to subtract the lake/reservoir storage from the GRACE time series (removal approach) and (b) to relocate the surface water storage at its exact location of origin (relocation approach). The earthquake correction (RECOG-EQ) was estimated from GRACE monthly solutions using the Bayesian approach provided in Einarsson et al. (2010) and takes into account both the co-seismic and the post-seismic signal.

This paper is organized as follows: Sect. 2 describes the input data and processing steps, presents the resulting correction products, and visualizes some of its key characteristics. In Sect. 3 we then show different exemplary applications of the dataset in order to illustrate its impact and to demonstrate its value: influence of RECOG on GRACE time series including a short discussion of the impact on data assimilation into a global hydrological model, the detection of drought indices in an earthquake-affected region, and a validation of the correction product using GNSS-observed station displacements. This is followed by a discussion of the benefits and limitations of the correction product in Sect. 4, before Sect. 5 summarizes the findings and gives an outlook on further development options.

2 Methods and results

In this section we describe the various input data and their sources (Sect. 2.1.1 and 2.1.2) that were used to perform the forward modelling (Sect. 2.1.4) of the surface waterbodies and its necessary pre-processing steps (Sect. 2.1.3). The resulting lake/reservoir correction dataset, RECOG-LR, is discussed in Sect. 2.1.5. The processing steps for computing the earthquake correction are described in Sect. 2.2.1, followed by the results of RECOG-EQ (Sect. 2.2.2).

2.1 Lake/reservoir correction, RECOG-LR

The lake/reservoir correction is based on a subset of currently 283 of the largest surface waterbodies monitored with satellite altimetry (see Supplement S1 for a complete list). The lake water volume variations product is designed around (1) monthly water level time series from a global multi-satellite product and (2) the (static) surface water extent area for each lake.

2.1.1 Lake level time series from altimetry

Satellite altimetry measures the distance between the satellite and the Earth surface (i.e. the range) by analysing the transmitted and received radar echo after it has been reflected by the Earth's surface. Originally, the technique was developed for open-ocean applications. However, if the data are carefully pre-processed, they can also be used for estimating the height of inland waterbodies, such as lakes and reservoirs. Since the inland signals are frequently contaminated by land reflections, a rigorous outlier detection (Schwatke et al., 2015) as well as a dedicated retracking (e.g. Gommenginger et al., 2011; Passaro et al., 2018) is mandatory.

In this study, time series created by DAHITI (Schwatke et al., 2015) are used. DAHITI provides water level time series of more than 2000 globally distributed inland targets (i.e. lakes, rivers, and reservoirs) in a period between 1992 and today, depending on the satellite mission covering the waterbody. Data from altimeter missions TOPEX, Jason-1/-2/-3, ERS-2, Envisat, SARAL, Sentinel-3A/-3B, ICESat, and Cryosat-2 are combined in a Kalman filter approach after being retracked with the Improved Threshold retracker (Bao et al., 2009). A key element of the approach is an extended outlier detection before and after data combination, including an optional classification of the radar echoes. The temporal resolution of the time series differs depending on the size of a lake. Small lakes that are only covered by one single satellite track can be measured every 35 or 10 d (depending on the mission), whereas for large lakes a height can be derived almost every day. Moreover, information can only be provided for those waterbodies located directly beneath a satellite's tracks (Dettmering et al., 2020), preventing the creation of water level time series of small lakes located between the satellites' ground tracks. In addition, for small lakes or lakes surrounded by large topography no reliable height information might be created due to corrupted or too-noisy radar echoes.

The quality of the DAHITI water level time series depends on various criteria, mainly on the size of the lake and the length of the crossing satellite track as well as the surrounding topography. Comparison with in situ data show RMSE of a few centimetres for larger lakes and RMSE of some decimetres for river crossings (Schwatke et al., 2015).

2.1.2 Creating lake shapes from remote sensing

Based on MODIS optical satellite data, daily surface water extents are provided by the DLR's Global WaterPack (GWP) product (Klein et al., 2017). To receive reliable estimates of the extent of large global lakes and reservoirs, daily observations are aggregated to obtain maximum waterbody extents for the years 2003 to 2018. To capture coherent waterbodies, a pixel-based region-growing algorithm is applied, using ancillary information of the temporal static HydroLAKES dataset (Messager et al., 2016) for waterbody identification. Hereby, every water pixel in the aggregated GWP raster layer that spatially overlaps with the original HydroLAKES shape file is assigned to the lake ID given by the HydroLAKES database. Subsequently, a seed point in every designated waterbody is determined, from which 8-pixel growth of the search window region is initiated, thus identifying neighbouring water pixels. This ensures that waterbodies are represented by coherent pixel groups only. After the growing process is finished, results are vectorized. With this dynamic approach, the risk of over- or underestimation of the actual water surface extent is reduced (see Fig. S2 in the Supplement for further details).

2.1.3 Data pre-processing

The input data were combined taking into account several pre-processing steps: water level observations were averaged to a monthly mean for consistency with the temporal resolution of the GRACE gravity field models. Missing months were linearly interpolated. The water level time series were cut to the investigation period (January 2003–December 2016) and then reduced by their respective means. To ensure a quick update of the correction product when new lakes will be added to the source databases or when their time series will be updated, the algorithm for matching water level time series with their respective lake surface area (as well as most of the following workflow) was automated. The first step for the combination was an automatically generated data table with global common lake IDs. Where no IDs were available, matching was achieved by comparing names while making sure that no double naming occurred in the input data. If that was not possible or no names were given, matches were found using a very strict search algorithm for the nearest lake shape to a given time series. If none of the above methods was successful, the respective lake was dismissed and not included in the correction. In total, matches for 283 lakes/reservoirs were found for RECOG-LR RL01 (see Fig. 2); a detailed list is provided in the Supplement (Sect. S1). The surface waterbody shapes were then discretized on a fine-resolution, 0.025^∘ grid to be able to capture long but narrow reservoirs in valleys. As in our algorithm grid cells are only assigned to belong to a waterbody if their midpoint lies within the lake polygon, small waterbodies would otherwise be missed completely. Multiplication with the altimetry-derived water height resulted in water volume estimates for each of the grid cells. These water volumes were subsequently distributed proportionally over a lower-resolution, 0.5^∘ grid to save computation time in the forward modelling and because a 0.5^∘ resolution is more than sufficient for applications to GRACE. The resulting global grids of lake/reservoir-related water height anomalies for each of the 168 months of our investigation period then entered the forward-modelling algorithm.

Figure 2 Overview of all 283 lakes/reservoirs in the dataset (blue areas) given on a 0.5^∘ grid.

2.1.4 Forward modelling

The localized altimetry/remote-sensing-derived surface water variations have to be converted to the GRACE spatial resolution before they can be subtracted from monthly GRACE gravity field estimates. In this forward-modelling step, the gridded values were expanded into spherical harmonic coefficients up to degree (n) and order (m) 96 according to

$$\begin{array}{ll}\text{(1)}& {\displaystyle }\left[\begin{array}{c}\mathrm{\Delta }{C}_{nm}\\ \mathrm{\Delta }{S}_{nm}\end{array}\right]=& {\displaystyle }{\displaystyle \frac{R^{\mathrm{2}}}{M}}\cdot {\displaystyle \frac{k_n+\mathrm{1}}{\mathrm{2}n+\mathrm{1}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}(\mathrm{\Delta }\text{TWS}\left(\mathit{\theta },\mathit{\lambda }\right){\displaystyle }& {\displaystyle }\cdot P_{nm}\left(\cos \mathit{\theta }\right)\cdot \left[\begin{array}{c}\cos \left(m\mathit{\lambda }\right)\\ \sin \left(m\mathit{\lambda }\right)\end{array}\right]\cdot \sin \left(\mathit{\theta }\right))\cdot \mathrm{d}\mathit{\lambda }\cdot \mathrm{d}\mathit{\theta },\end{array}$$

with ΔC_nm and ΔS_nm being the spherical harmonic coefficients (at degree n and order m), R the radius of the Earth, M the mass of the Earth, k_n the load Love numbers (Farrell, 1972), ΔTWS(θλ) the changes in altimetry-derived water storage in relation to colatitude θ and longitude λ, and P_nm the Legendre functions. Subsequently a standard GRACE filter was applied for smoothing (DDK3; Kusche, 2007; Kusche et al., 2009). The filtered coefficients are then denoted by $\mathrm{\Delta }{C}_{nm}^{\mathrm{F}}$ and $\mathrm{\Delta }{S}_{nm}^{\mathrm{F}}$. Differently filtered corrections are available upon request. This gives us the idealized signal that GRACE would measure if it were influenced by the changing mass in the lakes/reservoirs only. For a grid-based evaluation ($\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid) a recomputation using Eq. (2) is necessary to calculate the lake water storage ΔTWS^F(θ,λ) for every grid cell after filtering (Wahr et al., 1998), again up to degree and order 96:

$$\begin{array}{ll}\text{(2)}& {\displaystyle }\mathrm{\Delta }{\text{TWS}}^{\mathrm{F}}\left(\mathit{\theta },\mathit{\lambda }\right)=& {\displaystyle }{\displaystyle \frac{M}{\mathrm{4}\mathit{\pi }{R}^{\mathrm{2}}\mathit{\rho }}}\sum _{n=\mathrm{0}}^{\mathrm{96}}\sum _{m=\mathrm{0}}^n({\displaystyle \frac{\mathrm{2}n+\mathrm{1}}{\mathrm{1}+k_n}}\cdot P_{nm}(\cos \mathit{\theta }){\displaystyle }& {\displaystyle }\cdot \left(\mathrm{\Delta }{C}_{nm}^{\mathrm{F}}\cos \left(m\mathit{\lambda }\right)+\mathrm{\Delta }{S}_{nm}^{\mathrm{F}}\sin \left(m\mathit{\lambda }\right)\right)).\end{array}$$

Here we included the density of water ρ (1025 kg m⁻³) to obtain a lake water storage result in metres of equivalent water heights (EWHs), corresponding to the input variations in water storage from the water level time series.

2.1.5 Results for RECOG-LR

The resulting product for the lake and reservoir correction consists of two different parts, namely (1) the forward-modelled lake water correction to be subtracted from the GRACE data to remove the influence of lakes/reservoirs (removal approach). It is provided both on a spherical harmonic level (coefficients $\mathrm{\Delta }{C}_{nm}^{\mathrm{F}}$ and $\mathrm{\Delta }{S}_{nm}^{\mathrm{F}}$) and as a gridded data product (ΔTWS^F(θ,λ) from Eq. 2). The second part consists of (2) the altimetry-derived monthly water levels for each 0.5^∘ grid cell that can be used to re-add the measured lake volume to its actual area (relocation approach). Figure 2 highlights each grid cell that includes surface waterbodies with data used for the correction. Note that, although most of the major lakes and reservoirs are covered, some had to be excluded from the dataset for reasons explained further in Sect. 4. Figures showing an exemplary seasonal cycle of RECOG-LR can be found in the Supplement (Fig. S3), and an animation of the monthly changes of the lake/reservoir water storage for the full time series is provided in the “Video supplement” section (Deggim et al., 2020b; https://doi.org/10.5446/48188).

Figure 3 RECOG-LR on a global scale with its (a) mean seasonal amplitude and (b) trend.

Figure 3a shows the mean amplitudes of the seasonal variations of the correction for each grid cell. The most prominent features are the Caspian Sea in Asia and the Great Lakes in North America with amplitudes of about 10 cm. However, peaks for individual months can reach as much as 30 cm of TWS correction, as shown for exemplary time series in Fig. S3. Most of the other lakes and reservoirs have mean correction amplitudes in the area of 0 to 3 cm. Figure 3b displays the linear trend of the lake correction product. Again, the most distinctive features are a strong negative trend of the Caspian Sea and a strong water storage increase in the Great Lakes (altimetry time series shown in Fig. S4 in the Supplement for comparison). Strongly visible is also a positive trend (mass increase) in Lake Victoria accompanied by a clear mass loss in the nearby Lake Malawi. However, it also becomes evident that some surface waterbodies can have a rather prominent trend signal without having a strong seasonality, as for example Lake Oahe in South Dakota, or exhibit a strong seasonality without any long-term trend, such as Lake Chad in Africa or Lake Guri and Tucurui Reservoir in South America. In other surface waterbodies, such as the artificial reservoir Lake Volta, the time series is dominated by a strong inter-annual signal (Ni et al., 2017), which does not show up prominently in Fig. 3 but is very much visible in the total temporal variability shown in Fig. 5 below.

2.2 Earthquake correction, RECOG-EQ

2.2.1 Processing of RECOG-EQ

Different from the lake/reservoir correction, which is computed by forward modelling using independent datasets (altimetry/remote sensing), the earthquake correction is derived by fitting a parametric function to monthly GRACE data along the following processing line: in a first step, spherical harmonic coefficients have to be backward-modelled to gridded geoid changes by

$$\begin{array}{ll}\text{(3)}& {\displaystyle }\mathrm{\Delta }\text{GC}\left(\mathit{\theta },\mathit{\lambda }\right)=& {\displaystyle }R\sum _{n=\mathrm{0}}^{\mathrm{96}}\sum _{m=\mathrm{0}}^n\left(P_{nm}\right(\cos \mathit{\theta }){\displaystyle }& {\displaystyle }\cdot \left(\mathrm{\Delta }{C}_{nm}\cos \left(m\mathit{\lambda }\right)+\mathrm{\Delta }{S}_{nm}\sin \left(m\mathit{\lambda }\right)\right))\end{array}$$

to be able to apply the Bayesian approach provided in Einarsson et al. (2010). The total geoid changes for a specific location (θ_i,λ_j) can be subdivided into a bias (ΔGC_bias), trend (ΔGC_trend), annual- (ΔGC_ann) and semi-annual signal (ΔGC_semiann), S2 aliasing period of 161 d (ΔGC_N2), and the earthquake signal (ΔGC_EQ) as

$$\begin{array}{ll}\text{(4)}& {\displaystyle }\mathrm{\Delta }\text{GC}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j,t\right)=& {\displaystyle }\mathrm{\Delta }{\text{GC}}_{\text{bias}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j,t\right){\displaystyle }& {\displaystyle }+\mathrm{\Delta }{\text{GC}}_{\text{trend}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j,t\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{\Delta }{\text{GC}}_{\text{ann}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j,t\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{\Delta }{\text{GC}}_{\text{semiann}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j,t\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{\Delta }{\text{GC}}_{\text{S2}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j,t\right)+\mathrm{\Delta }{\text{GC}}_{\text{EQ}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j,t\right),\end{array}$$

which contains the model coefficients C_bias, C_trend, C_ann, ϕ_ann, C_semiann, ϕ_semiann, C_S2, and ϕ_S2. The earthquake signal included here is described by a co-seismic and a post-seismic component modelled as

$$\begin{array}{ll}\text{(5)}& {\displaystyle }\mathrm{\Delta }{\text{GC}}_{\text{EQ}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j,t\right)=& {\displaystyle }{C}_{v_{\text{co}}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j\right)H_{t_v}\left(t\right){\displaystyle }& {\displaystyle }+C_{v_{\text{post}}}\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j\right)H_{t_v}\left(t\right)\\ {\displaystyle }& {\displaystyle }\cdot \left(\mathrm{1}-e^{-\frac{t-t_v\left({\mathit{\theta }}_i,{\mathit{\lambda }}_j\right)}{\mathit{\tau }({\mathit{\theta }}_i,{\mathit{\lambda }}_j)}}\right).\end{array}$$

$C_{v_{\text{co}}}$ and $C_{v_{\text{post}}}$ describe coefficients for the co- and post-seismic component of the respective earthquake v, $H_{t_v}\left(t\right)$ is the Heaviside step function at time t_v, and τ is the decay rate. All coefficients are then estimated using Monte Carlo integration for quasi-linear models and are used to estimate the total earthquake signal. This signal is then removed from the total geoid changes for each considered earthquake to derive an earthquake-corrected dataset. Furthermore, we applied a spatial radial Gaussian window with a radius (half width) of 157 km to consider only regions that were affected by earthquakes. The centre of the Gaussian window is placed in the epicentre of the respective earthquake. Einarsson's approach is, as recommended, consecutively applied to earthquakes with a magnitude that is larger than or equal to 9.0, which is a criterion met by the Sumatra–Andaman earthquake (M9.1) in December 2004 and the Tohoku earthquake (M9) in March 2011. Einarsson (2011), for example, showed that earthquakes with a magnitude larger than 9.0 are clearly visible in GRACE data, while earthquakes with lower magnitude cannot always be clearly separated. For more information about the approach see Einarsson et al. (2010) and Einarsson (2011). To derive TWS anomalies, the geoid changes are forward-modelled to spherical harmonic coefficients similar to Eq. (1) by

$$\begin{array}{ll}\text{(6)}& {\displaystyle }\left[\begin{array}{c}\mathrm{\Delta }{C}_{nm}\\ \mathrm{\Delta }{S}_{nm}\end{array}\right]=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{R}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}(\mathrm{\Delta }\text{TWS}\left(\mathit{\theta },\mathit{\lambda }\right)\cdot P_{nm}(\cos \mathit{\theta }){\displaystyle }& {\displaystyle }\cdot \left[\begin{array}{c}\cos \left(m\mathit{\lambda }\right)\\ \sin \left(m\mathit{\lambda }\right)\end{array}\right]\cdot \sin \left(\mathit{\theta }\right))\cdot \mathrm{d}\mathit{\lambda }\cdot \mathrm{d}\mathit{\theta }\end{array}$$

and again backward-modelled to TWS changes as described in Eq. (2). The final earthquake correction product, RECOG-EQ, is then derived by computing the difference between the uncorrected and earthquake-corrected TWS anomalies (TWSA).

2.2.2 Results for RECOG-EQ

The correction is provided equivalently to the lake correction: the dataset is processed on a global 0.5^∘ grid using spherical harmonic coefficients up to degree and order 96, is DDK3-filtered (different filters are available upon request), and covers the period 2003 to 2016.

Figure 4 Linear trends (January 2003–December 2016) of TWS changes (cm EWH yr⁻¹) of the earthquake correction that includes the (a) Sumatra–Andaman earthquake from 2004 and the (b) Tohoku earthquake from 2011.

The correction only shows differences over the regions of the 2004 Sumatra–Andaman earthquake (Fig. 4a) and the 2011 Tohoku earthquake (Fig. 4b) because we only corrected for these two earthquakes. For the Sumatra–Andaman region the linear trends of the correction reach from about −2.7 to 1.1 cm EWH per year. Negative linear trends of down to −2.7 cm EWH per year can be found north of Indonesia and west of peninsular Malaysia, while positive trends are visible in the Indian Ocean close the coast of Sumatra. Considering Tohoku, the linear trends range from about −2.5 to 1.4. The dominant negative part can be found to the west of the Tohoku region, while the positive parts are apparent in the Pacific Ocean, southeast of Tohoku. These results suggest that uncorrected TWS changes might hinder the correct analysis of the data for hydrological studies, because the post-seismic part of the earthquake might falsely be interpreted as a linear trend in the uncorrected TWS changes. This is particularly relevant when the earthquake occurs at the beginning of the time series as is the case for the Sumatra–Andaman earthquake. The results shown here are derived from the ITSG-Grace2018 solutions (Kvas et al., 2019). Earthquake corrections derived from other GRACE solutions provide similar findings; for completeness they are attached in the Supplement (Figs. S5 and S6). We do not apply geophysical forward modelling since these models heavily depend on dislocation parameters, fault geometry, and background rheology, and parameters are typically tuned to fit seismic and GNSS measurements but do not fit observed gravity changes.

3 Applications and validation

In this section we would like to show the influence and benefit of the correction datasets. For this purpose, in Sect. 3.1 we first illustrate the impact of subtracting the two corrections RECOG-LR and RECOG-EQ from the GRACE time series in terms of change in signal variability and trends. This includes a short discussion of the benefit of RECOG-LR for data assimilation into hydrological models, one of the main target applications of the corrected GRACE dataset. The comparison to GNSS surface displacements (Sect. 3.2) not only shows the influence of the surface water correction on geometrical surface deformations even several hundreds of kilometres around lakes/reservoirs but also represents a first validation of the corrected GRACE time series using independent observations. Finally, the influence of removing earthquake signals (RECOG-EQ) for hydrological drought detection is described in Sect. 3.3.

3.1 Influence of RECOG RL01 on a global GRACE time series

3.1.1 Influence of lake/reservoir correction, RECOG-LR, on GRACE

We first investigate the influence of subtracting the lake/reservoir correction dataset from a global GRACE time series. For this purpose, we derive gridded TWSA from the ITSG-Grace2018 (Mayer-Gürr et al., 2018; Kvas et al., 2019) spherical harmonic expansion up to degree and order 96 considering corrections for low-degree coefficients (Swenson et al., 2008; Cheng et al., 2011; Sun et al., 2016), glacial isostatic adjustment (A et al., 2013), and applying the DDK3 filter (Kusche, 2007).

For the removal approach, we then reduce the lake correction (Sect. 2.1.5) from the GRACE-derived TWSA grids. For the relocation approach, we re-add the altimetry-derived monthly water mass estimates. This leads to three global TWSA datasets: (1) GRACE-based only, (2) GRACE TWSA after removing altimetry-based lake/reservoir storage, and (3) GRACE TWS with relocated altimetry-based lake signal.

Figure 5 Temporal root mean square (rms) of the lake/reservoir correction time series for each grid cell (a) and the relative reduction in temporal rms when subtracting the correction (removal) from the original GRACE TWS time series (b).

Figure 5a shows the temporal root means square (rms) of the lake correction for each grid cell. Subtracting this correction reduces the temporal rms variability in the GRACE time series (Fig. 5b) by up to 75 % around the Caspian Sea in Asia and 50 % around Lake Victoria in Africa, compared to the original variability in ITSG-Grace2018. Values around the other lakes vary between 0 and 30 % with a few negative values in the area south of the Caspian Sea and in Canada. The latter can most likely be attributed to Gibbs oscillations by the bandlimited spectral representation of the data and/or to being spuriously introduced by the non-isotropic DDK3 filter applied to the lake signals. However, the option that the lake signal was hiding another impact, which can only be recovered after subtracting the correction, should not be completely ruled out, e.g. in the case of water transfer between compartments as from glacier/snow to lake water (Castellazzi et al., 2019).

Figure 6 Linear trend of GRACE TWS anomalies (January 2003–December 2016) in the Caspian Sea region (a–c) and in the Mississippi Basin and around the Great Lakes (d–f) without any correction (a, d), after removing the influence of lakes (b, e), and after relocating the lake signal (c, f). Please note that the Aral Sea has been excluded from RECOG-LR RL01 due to its strongly varying surface area, which is not yet captured in the database.

Figure 6 shows the influence of the lake correction on the linear trend in the GRACE time series for two detailed examples; a global map can additionally be found in the Supplement (Fig. S7). In the area around the Caspian Sea (Fig. 6a–c), a very strong negative trend of around −3 cm yr⁻¹ in the original GRACE time series (a) is almost completely removed by the lake correction (b). The relocation approach then restores the altimetry-derived lake water variation to the lake area (c). The second example shows the Mississippi Basin (Fig. 6d–f). Even though the Great Lakes are not part of the basin, they still have an effect, particularly on subbasin Alton (upstream from Alton, Illinois, USA), which is closest to the Great Lakes. Influences of the lake variations on the GRACE basin average of this subbasin can reach values of up to 5 cm in TWS for some months. This example also shows that a positive trend visible in the original GRACE data can mainly be attributed to surface water change and can be levelled by the correction. Subtracting the lake correction also reduces the positive trend partly caused by smaller lakes/reservoirs (e.g. Lake Oahe and Lake Sakakawea) in the Hermann subbasin in the northwest of the Mississippi.

The importance of correcting the signal in the nearby Great Lakes and the smaller lakes/reservoirs within the subbasins can also be detected when assimilating GRACE-derived TWSA into a hydrological model. The original and the RECOG-corrected, subbasin-averaged GRACE TWSA of the Mississippi River basin were assimilated into the Water – Global Assessment and Prognosis (WaterGAP; Döll et al., 1999; Alcamo et al., 2003; Müller Schmied et al., 2014) hydrological model, which has a $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid spatial resolution. Our assimilation theory follows Eicker et al. (2014) and Schumacher et al. (2015, 2016); for more information see the Supplement (Sect. S8).

In the Alton subbasin, the assimilation of the original GRACE observations introduces a spurious positive mass trend of 2.5 mm EWH yr⁻¹ not present in the original WaterGAP version. This trend is assumed to not originate from storage increase within the subbasin itself but to be caused by leakage due to a water storage increase in the Great Lakes (see Fig. 6), particularly in the nearby Lake Superior and Lake Michigan. Applying the RECOG correction dataset before assimilation reduces the linear trend in the subbasin to 0.9 mm EWH yr⁻¹ and thus prevents the spurious trend from appearing in the model output. This can also be confirmed on the grid cell level with a trend reduction in 53 % of the grid cells in the Alton subbasin, 40 % of them by more than 80 % compared with the results of assimilating the original TWSA. The strongest reductions appear in the northeastern part of the basin, e.g. from 12.3 to 3.5 mm EWH yr⁻¹ in a grid cell in close proximity to the Great Lakes (lat 46.25^∘, long −90.25^∘). In grid cells directly affected by lake water storage the effect can be even larger; e.g. for Lake Sakakawea (lat 48.25^∘, long $-\mathrm{103.25}{}^{\circ }$) the trend in the assimilation results changes from 493.5 to 8 mm EWH yr⁻¹, and the signal rms from 2548 to 120 mm. The results reveal the importance of applying a lake water correction to GRACE data before data assimilation not only for areas covered by waterbodies but also for neighbouring grid cells affected by strong leakage-in effects.

3.1.2 Influence of earthquake correction, RECOG-EQ, on GRACE

This section presents the application of the earthquake correction (Sect. 2.2) to the GRACE data. Linear trends for the period 2003 to 2016 are derived from (1) the original GRACE TWSA and (2) the corrected TWSA after applying RECOG-EQ. The correction changes the spatial pattern of the linear trend in the Sumatra–Andaman region (Fig. 7a and c), and the magnitude of positive trends increases by about 0.5 cm EWH yr⁻¹ from 4.2 to 4.7 cm EWH yr⁻¹. The spatial extension of the positive trends in the corrected data (Fig. 7c) reaches to peninsular Malaysia. Thus, the former slightly negative trends of about −1 cm EWH yr⁻¹ identify a smaller magnitude or even slightly positive values in the earthquake-corrected dataset. The change in trends might bias the correct analysis of linear trends for this region.

Figure 7 Linear trends (January 2003–December 2016) of TWS changes (cm EWH yr⁻¹) before and after removing earthquake signals of the Sumatra–Andaman earthquake from 2004 (a, before; c, after) and the Tohoku earthquake from 2011 (b, before; d, after).

When analysing the Tohoku earthquake results (Fig. 7b and 7d), we also see that the magnitude and spatial pattern of the trends change. In this case, the difference between original and corrected GRACE data is more obvious than in the Sumatra–Andaman region: the original dataset shows positive linear trends of about 2 cm EWH yr⁻¹ in the west of Tohoku, while negative trends of about −2 cm EWH yr⁻¹ can be found in the Pacific Ocean along the southeastern coast. These trend signals vanish almost completely after subtracting the earthquake correction, confirming our assumption made in Sect. 2.2.2: the post-seismic earthquake component was identified as linear trends in the original GRACE data and has clearly biased the trend analysis, leading to misinterpretation of the trends, especially for the Tohoku region.

Figure 8 Spatially averaged GRACE TWS anomalies for the original (black) and earthquake-corrected (blue) GRACE data and its correction (red) for (a) peninsular Malaysia (West Malaysia) and (b) Japan.

To analyse the effect of earthquakes on TWS changes in more detail and over land, spatially averaged TWS anomalies are compared for peninsular Malaysia in Fig. 8a and for Japan in Fig. 8b. Additional figures showing the signal in the epicentre of the two earthquakes are shown in the Supplement (Fig. S9). Regarding peninsular Malaysia, the uncorrected TWSA (black) show a strong decrease in TWS changes beginning in 2004, which results from the Sumatra–Andaman earthquake. After applying the correction (blue), this strong decrease is removed. The correction (red) shows nicely the co-seismic component of about −2.5 cm EWH as a jump between December 2004 and January 2005 and a following post-seismic relaxation, which increases the total correction towards −6 cm. Similar findings can be observed for the results for Japan and the Tohoku earthquake: the correction contains a co-seismic component down to −6 cm with post-seismic relaxation afterwards, but in this case the relaxation is slowly decreasing the amount of correction. However, both examples clarify that uncorrected GRACE data could lead to wrong conclusions about the underlying signals as the apparent trends can, for instance, hamper the identification of drier and wetter years in the GRACE time series.

3.2 Validation of RECOG-LR with GNSS

Global Navigation Satellite System time series contain signals resulting from surface mass redistribution: atmospheric pressure, non-tidal oceanic loading, and hydrological changes. Vertical displacement due to hydrological loading can be predicted by calculating the elastic response of an Earth model to the TWS changes. Previous studies have shown good agreement between modelled deformation from GRACE TWS and vertical displacement observed by GNSS (e.g. Springer et al., 2019; Tregoning et al., 2009; van Dam et al., 2007), while horizontal displacements due to hydrological loadings are much smaller than the vertical ones and have been found to be systematically underpredicted and out of phase (Chanard et al., 2018).

Here we use vertical displacement residuals from the ITRF2014 stacking (Altamimi et al., 2016) resulting from the second reprocessing campaign (repro2) by the International GNSS Service (IGS) (Rebischung et al., 2016) to validate the GRACE lake/reservoir correction data product RECOG-LR. The station velocities and discontinuities have been carefully removed from the GNSS time series. To be consistent with GRACE TWS, the effects of atmospheric and non-tidal oceanic loading are subtracted from the time series using the AOD1B product (Dobslaw et al., 2017). The displacements due to non-tidal oceanic and atmospheric loading are calculated at daily epochs to be consistent with the GNSS observations. GRACE TWS anomalies are converted into displacements in the spatial domain using spatial convolution of a point mass load Green's function (Farrell, 1972) in the centre-of-figure frame using a set of high-degree load Love numbers determined by H. Wang et al. (2012) for the Preliminary Reference Earth Model (PREM) (Dziewonski and Anderson, 1981). Finally, the daily GNSS displacements are averaged to monthly intervals before removing the temporal mean and linear trend from both GNSS-observed and GRACE-derived (modelled) displacements. Exemplary time series for GNSS-observed and GRACE-modelled displacements at different stations in the Great Lakes, Lake Victoria, and Caspian Sea regions are shown in the Supplement (Fig. S10). Here the GNSS time series are used to assess the impact of the lake/reservoir correction on GRACE data. To get an idea of the magnitude of its influence, we first show the temporal rms of the forward-modelled RECOG-LR signal (Fig. 9a) for ITRF2014 GNSS sites around the Great Lakes. This represents the signal variability of modelled station displacements caused only by the surface water variations described in RECOG-LR. The rms of the vertical displacements amounts to up to 1.4 mm for station ALGO (∼190 km away from the lake shore) and can reach higher values for other stations provided by NGL (Blewitt et al., 2018) (not shown), such as 2.3 mm at station BAKU close to the Caspian Sea and 2.5 mm at CHB5 directly at the shore of Lake Huron.

Figure 9 Temporal rms of the vertical displacement caused by the forward-modelled lake/reservoir signal of RECOG-LR (a), reduction of rms of GNSS time series after subtracting modelled deformation (b), and difference in rms reduction when subtracting GRACE after removing and relocating RECOG-LR vs. using the uncorrected GRACE signal (c).

To investigate the agreement of observed (GNSS) and GRACE-derived displacement time series, we compute the reduction in temporal rms after subtracting the modelled displacements from the observed time series. For this we first compute the signal variability rms_GNSS of the observed station displacements and subsequently the variability rms_GNSS-GRACE of the GNSS time series after subtracting GRACE-derived station movements. The relative reduction in the rms is given by the following equation:

$$\begin{array}{}\text{(7)}& {\text{rms}}_{\text{reduction}}\left(\mathit{\%}\right)={\displaystyle \frac{{\text{rms}}_{\text{GNSS}}-{\text{rms}}_{\text{GNSS-GRACE}}}{{\text{rms}}_{\text{GNSS}}}}\times \mathrm{100}.\end{array}$$

This relative reduction in rms is shown in Fig. 9b, and it can be observed that, by removing the modelled displacement based on GRACE from the GNSS time series, all stations around the Great Lakes have their rms reduced, indicated by positive rms reductions. This means for ITRF2014 GNSS sites around the Great Lakes we can explain some of the observed vertical displacement by TWS changes at most of the stations, with a stronger agreement at stations closer to the Great Lakes (ALGO and NRC1 stations).

To assess the influence of RECOG-LR on these reductions, we compute two different versions of Eq. (7): we compare the rms reductions using (1) the original (unmodified) GRACE data and (2) GRACE after subtracting and relocating the lake signal. The difference between these two rms reductions is plotted in Fig. 9c. A positive value means that the RECOG-corrected GRACE signal explains a better portion of the observed GNSS station displacements than the uncorrected GRACE signal. All but one of the stations around the Great Lakes show a positive effect of the lake correction. The only station unimproved is located quite far from the lakes, directly at the coast of the ocean, and might primarily be influenced by oceanic leakage as discussed in van Dam, et al. (2007). The largest improvements can be observed at stations ALGO (6.1 %) and NRC1 (4.3 %), both around 180–200 km away from Lake Huron. In the time series plot in Fig. S10 of the Supplement this improvement is indicated by a better fit of the corrected GRACE time series with the GNSS observations. Also for other ITRF2014 stations in the area of large surface waterbodies we find a systematic improvement of the rms reduction when applying RECOG-LR. Examples are stations around the Caspian Sea (NSSP: 0.2 % improvement; TEHN: 0.3 % improvement) and around Lake Victoria (NURK: 1.4 % improvement; RCMN: 2.7 % improvement). Here it should be noted that these stations are not located close to the shore but between 100–400 km away from the lakes, yet they are still affected by the correction of the lake signal. To put these numbers of up to 6 % improvement into perspective, a change in the Earth model used for the conversion of TWS to deformation, which has previously been found to be relevant for GNSS analysis (Karegar et al., 2017), has a much smaller influence. The differences in rms reduction using, for example, different sets of load Love numbers amount to only <0.4 %. Thus, from the above results, we are confident that correcting for the leakage effect of lake/reservoir water storage in GRACE time series can have a considerable positive effect on the comparison of GRACE and GNSS observations.

3.3 Hydrological drought detection with earthquake correction, RECOG-EQ

Hydrological drought detection using GRACE data has been applied in various studies, for example in Houborg et al. (2012), Thomas et al. (2014), Zhao et al. (2017), Boergens et al. (2020), and Gerdener et al. (2020a) for many regions of the Earth. Removing the earthquake signal before the analysis might lead to a change in the drought detection results, which could have a significant impact on, for example, the decisions of policy makers. In this section, we show an example of the influence of the earthquake correction (Sect. 2.2) for detecting drought events in peninsular Malaysia. To show the longer-term behaviour of droughts, the drought severity index using accumulation (DSIA) used in Gerdener et al. (2020a) is computed. As typically done with meteorological indicators, the observable is accumulated for a chosen period q ($\mathrm{\Delta }{\text{TWS}}_{i,j,q}^{+}$) before its computation because we refer it to a duration of drought. For example, if the accumulation period q is set to 3 months, the accumulated TWS changes for March 2003 will be the sum of TWS changes of January, February, and March. The accumulation also serves as a running mean. The DSIA is then computed by

$$\begin{array}{}\text{(8)}& {\text{DSIA}}_{i,j,q}={\displaystyle \frac{\mathrm{\Delta }{\text{TWS}}_{i,j,q}^{+}-\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\text{TWS}}_{j,q}^{+}}}{{\mathit{\sigma }}_{j,q}^{+}}},\end{array}$$

where $\mathrm{\Delta }{\text{TWS}}_{i,j,q}^{+}$ is the accumulated TWS changes in year i and month $j=\mathrm{1},\mathrm{\dots },\mathrm{12}$; $\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\text{TWS}}_{j,q}^{+}}$ is the mean monthly accumulated TWS change, e.g. the mean over all Januaries; and ${\mathit{\sigma }}_{j,q}^{+}$ is the monthly standard deviation. Here, it is used to identify hydrological drought events for an accumulation period of 6 months (DSIA6) over peninsular Malaysia. Figure 10 shows the resulting DSIA6 time series using corrected (blue) and uncorrected (black) TWS changes. The uncorrected DSIA6 identifies a mainly moderate dry period in 2010 and 2011 and a severe period at the beginning of 2014. Both periods are also identified using the corrected DSIA6 but with different intensity. The period 2010/11 is slightly more intense now, and the drought in 2014 is extreme (e.g. Tan et al., 2017). Furthermore, the corrected DSIA6 shows exceptional drought in 2005, which was not identified with the uncorrected DSIA6. These findings are supported by, for example, the EM-DAT database (EM-DAT, 2020) and Hashim et al. (2016), who also identified a drought over peninsular Malaysia in 2005. However, in March 2005, a second earthquake (Nias earthquake) also occurred close to the Sumatra–Andaman region with a magnitude of M8.6. As stated by, for example, Broerse (2014), earthquakes with a lower magnitude are not always clearly visible in the data, but it should be noticed that the Nias earthquake might still have possible influences on the time series analysis.

Figure 10 GRACE-derived TWSA DSIA6 over peninsular Malaysia (West Malaysia) shown without (black) and with (blue) earthquake correction.

4 Limitations

As mentioned in Sect. 2.1, static lake shapes were used to determine the lake volume and subsequently the mass estimates in this first release of RECOG. In contrast to monthly or daily lake shapes, these are not (or much less) impacted by cloud coverage in the optical satellite images and, thus, more reliable. This approximation works for most of the lakes with shores that are not too flat. Estimated errors caused by ignoring the changing surface area revealed numbers between 3–15 % for exemplary surface waterbodies (Semmelroth, 2019). Extreme cases with either very flat shores and high water level variations or strong trends such as the Aral Sea in Asia, whose surface area shrunk to a fraction of its original size in the last couple of decades, have been omitted in the correction.

Though RECOG-LR covers most of the major lakes around the world, some are not yet implemented in the dataset due to failed automatic matching between water level time series and lake surface area (e.g. Lake Athabasca, see also Sect. 2.1.3), insufficient water level time series due to ice coverage for major parts of the year (e.g. Lake Taymyr) or missing flyovers by altimeter satellites (e.g. Dead Sea; see also Sect. 2.1.1). More generally, we only consider surface waterbodies that can be captured by satellite altimetry and this way underestimate the impact of surface water storage in regions with a large number of small lakes and dams, e.g. in the US or in India.

Additionally, changes in surface water volume also affect surrounding groundwater storage (and thus GRACE data) by groundwater–surface water interactions (e.g. Bierkens and Wada 2019), which have not yet been considered in our data product. Steric effects caused by the thermal expansion of water have an influence on the altimetry-derived water levels in large surface waterbodies, which have not been accounted for in the conversion from water volume to mass change. Studies only exist for very large lakes, such as the Caspian Sea, where the steric water level change was found to amount to about $\mathrm{1}/\mathrm{3}$ of the seasonal amplitude (Chen et al., 2017; Loomis and Luthcke, 2017). We expect this to be an extreme case, with the influence on smaller waterbodies being considerably smaller. However, since the necessary temperature and salinity profiles are not available on a global scale including smaller waterbodies, no steric correction has been included in RECOG-LR.

5 Code availability

The GROOPS software toolkit has been used. It is available on GitHub: https://github.com/groops-devs/groops/commit/a52631fc1817acdc4b40e1caae546254f36a2653 (last access: 6 November 2020; Mayer-Gürr et al., 2020).

6 Data availability

RECOG includes (1) RECOG-LR with (1a) global gridded time series for the given time span with the correction values in total water storage for each grid cell and month, (1b) the same values in SH coefficients of degree and order 96, and (1c) the monthly gridded altimetry-derived water height for the relocation approach and (2) RECOG-EQ with the monthly gridded earthquake correction for TWS changes, containing the 2004 Sumatra–Andaman and the 2011 Tohoku earthquake. The data have been uploaded to the PANGAEA database: RECOG-LR: Deggim et al. (2020a), https://doi.org/10.1594/PANGAEA.921851; RECOG-EQ: Gerdener et al. (2020b), https://doi.org/10.1594/PANGAEA.921923.

7 Conclusions and outlook

Leakage effects of surface waterbodies and non-hydrology-related mass change signals have a strong influence on water storage estimates from GRACE, complicating its use for hydrological studies and specifically for calibration and data assimilation. Volume change estimates from combining satellite altimetry with remote sensing information can be used to remove (or at least strongly reduce) the effect of lakes and reservoirs from the GRACE data on a global scale, with particular benefit in regions close to big lakes/reservoirs or in regions with many smaller lakes and reservoirs. The earthquake signal (co-seismic and post-seismic), which masks hydrological variations in the vicinity of large earthquakes, can be extracted directly from the GRACE time series.

In this contribution, we introduced the first release of a new global correction dataset, RECOG RL01, for removing both the lake/reservoir storage (RECOG-LR) and the earthquake signal (RECOG-EQ) from the GRACE time series, while also offering the possibility to relocate the altimetry-derived mass change to its original surface waterbody outline.

Exemplary applications show that the correction product can reduce the signal variability (rms) of the GRACE signal by up to 75 % for the most prominent example of the Caspian Sea and that affected areas not only include the lake areas themselves but can also extend for tens to hundreds of kilometres around the waterbodies due to leakage. Special precaution has to be taken when assimilating GRACE data into hydrological models in the proximity of large surface waterbodies. In this context, the correction product is particularly valuable for models that do not include a surface water compartment at all, but the reduction of the leakage effect can also make it beneficial for models that do. For the example of the Alton subbasin of the Mississippi the leakage signal of the Great Lakes would cause an artificial mass increase in the assimilated model runs of WaterGAP, which can be prevented by subtracting and relocating the surface water storage before assimilation. A validation of the corrected GRACE signal using observed vertical GNSS station displacements shows an improvement of the fit between GRACE and GNSS of up to 6 % for stations at around 180–200 km distance from the Great Lakes. Applying the earthquake correction allows for the determination of several severe and one exceptionally severe drought events in peninsular Malaysia, which a GRACE-based drought indicator would otherwise have missed without first correcting for the earthquake signal. Therefore, depending on the application at hand, we recommend applying both RECOG-LR (globally) and RECOG-EQ (especially when research is performed in areas that underlie large earthquakes) as a standard post-processing step for analysing GRACE data.

Future improvements of the correction data product will introduce dynamic lake shapes replacing the static lake surface areas used so far, which will enable the inclusion of surface waterbodies with major area changes, such as the Aral Sea. Alternatively, pre-processed volume change time series of DAHITI (Schwatke et al., 2020) can be introduced for lakes for which they are already provided. The GRACE time series continues thanks to GRACE-Follow-On, and it is thus important to also provide a continuously updated surface waterbody correction which will rely on continuously updated source data. An extension of the RECOG-LR to include further surface waterbodies and a closure of existing data gaps as soon as new data are added to DAHITI will easily be possible thanks to a largely automated process of matching lake IDs in the altimetry database with corresponding lake shapes and forward-modelling the water volume to filtered GRACE-like TWS. Thus RECOG-LR can be updated as soon as new data become available.

So far, RECOG-LR focusses on lakes and reservoirs. However, there are other forms of surface waterbodies whose effects on GRACE data are still disregarded and not yet covered by any correction. An example for this are rivers, especially river deltas of big river basins (e.g. Mississippi, Amazon, Congo) that are highly influenced by strong seasonal variations in water flux as well as influences by tides in the estuary. Such a correction would be especially interesting for hydrological modelling.

The correction data product can also be extended to cover additional geophysical phenomena. For example, since most hydrological models do not include an explicit glacier compartment, it would be beneficial to extend the correction dataset to remove the glacier mass component from an independent glacier model or from remote sensing information before assimilating GRACE data into the model.

The presented correction products offer possibilities for more sophisticated data assimilation strategies. For example, GRACE data after removal of the lake/reservoir/earthquake correction can be assimilated into non-surface water compartments of a hydrological model, while the altimetry-derived relocation dataset can be assimilated solely into the surface water compartment of models that explicitly incorporate this. The best way to use this information for data assimilation will have to be further investigated.

Further research for the earthquake correction should consider comparing different methods and should also analyse possible influences of earthquakes with a magnitude lower than 9.0 in more detail.