Contributors: E. Boergens, C. Dahle, J. Haas, H. Dobslaw, F. Flechtner, A. Güntner

GFZ Helmholtz Centre for Geosciences

Issued by: GFZ / E. Boergens 

Date: 16/07/2025

Ref: C3S2_313c_EODC_WP3-DDP-TWSA-v1_202507_ATBD_i4

Official reference number service contract: ECMWF/COPERNICUS/2024/C3S2_313c_EODC

Table of Contents

History of modifications

Version 

Issue

Date 

Description of modification

Chapters / Sections

1.0

1

21/03/2025

Initial submission

All

1.0

2

23/04/2025

Document amended in response to independent review and finalised for publication.

All

1.0

3

12/06/2025

Small updates in Cost-G (Grace) table and Section 3.2.4, added Leakage to the general definitions

All

1.0

4

16/07/2025

Final links to COST-G RL02 data provided

Update of Figure 4 for COST-G RL02

Amending text in Section 2.1 for COST-G RL02

Small update in COST-G (GRACE) table in Section 2.1.1

Added links to PUGS

Changes in references and acronyms according to text changes

All

List of datasets covered by this document

Deliverable ID

Product title

Product type (CDR, ICDR)

Version number

Delivery date

WP3-CDR-TWSA-1

ECV Data Product  1

CDR

v1.0

30.06.2025

Acronyms

Acronym

Definition

AC

Analysis Centre

AIUB

Astronomical Institute of the University of Bern

ATBD

Algorithm Theoretical Basis Document

CCI

Climate Change Initiative

CF

Centre of Figure

CM

Centre of Mass

COST-G

Combination Service for Time-variable Gravity Fields

CSR

Centre for Space Research - University of Texas at Austin

C3S

Copernicus Climate Change Service

DDK

Proper name of a decorrelation and smoothing method for GRACE gravity fields, not an acronym

DLR

Deutsches Zentrum für Luft- und Raumfahrt (German Aerospace Center)

ESA

European Space Agency

GFZ

Helmholtz-Zentrum für Geoforschung (Helmholtz Centre for Geosciences)

GIA

Glacial Isostatic Adjustment

GPS

Global Positioning System

GRACE

Gravity Recovery and Climate Experiment

GRACE-FO

GRACE Follow-On

GSHHG

Global Self-consistent, Hierarchical, High-resolution Geography Database

HUST

HuaZhong University of Science and Technology

IAG

International Association of Geodesy

ICGEM

International Centre for Global Earth Models

IERS

International Earth Rotation and Reference Systems Service

IGFS

International Gravity Field Service

ILRS

International Laser Ranging Service

JPL

Jet Propulsion Laboratory - California Institute of Technology

KBR

K/Ka-Band Ranging

LAGEOS

Laser Geodynamics Satellite

LARES

Laser Relativity Satellite

NASA

National Aeronautics and Space Administration

PUGS

Product User Guide and Specification

SAR

Synthetic Aperture Radar

SDS

Science Data System

SH

Spherical Harmonic

SLR

Satellite Laser Ranging

Starlette

Satellite de Taille Adaptée avec Réflecteurs Laser pour les Études de la Terre

TUG

Graz University of Technology

TWS

Terrestrial Water Storage

TWSA

Terrestrial Water Storage Anomaly

VDK

Time-Variable DDK (see above)

General Definitions

Gravity Recovery and Climate Experiment (GRACE): Twin satellite mission observing the time variable gravity field of the Earth from which water mass redistribution can be inferred. The mission was active between 2002 and 2017 and was a joint US-German project. 

GRACE Follow-On (GRACE-FO): Identically constructed successor satellite mission of GRACE, launched in 2018. The mission is again a joint US-German project.

Terrestrial Water Storage Anomaly (TWSA): TWSA represents the deviation or anomaly of terrestrial water storage (TWS) in a certain epoch (e.g., a month) from the long-term mean TWS in all hydrological water storage compartments: groundwater storage, soil moisture, surface water storage, and snow storage and glaciers. Globally, TWSA can only be observed with GRACE and GRACE-FO. 

Spherical Harmonics: Mathematical functions defined on the surface of a sphere, commonly used to express potential fields such as the Earth's gravity field by solving Laplace's equation.

Spherical Harmonic Coefficients: Commonly used representation of a global gravity field model based on the Spherical Harmonics basis functions. They are jointly estimated by least-squares adjustment together with other (instrument and orbit) parameters using satellite observations as provided by GRACE and GRACE-FO.

Equivalent Water Height: Equivalent Water Height represents the mass change per unit area as if it were entirely due to a uniform layer of water. It is derived from changes in gravity and surface deformation and is typically expressed in millimetres or centimetres.

Uncertainty: Uncertainty measures jointly the precision and accuracy of the TWSA data, expressed as a standard deviation. The uncertainties of the TWSA grids are influenced, among others, by the uncertainties of the sensors on the satellites, uncertainties of background models applied in the processing chain, the orbit configuration of the missions (near-polar orbit), or environmental effects such as solar activity.

Glacial Isostatic AdjustmentGlacial isostatic adjustment describes the deformational behaviour of the solid earth in response to glacial loading processes, in particular resulting in surface displacements as in stress, gravity and sea level changes.

Satellite Laser Ranging: A space-geodetic technique measuring the run-time of a laser pulse from a ground station to a passive retroreflector on a satellite back to the ground station and thus allowing to determine the distance between the ground station and the satellite with millimeter level precision. SLR is mainly used to realize global terrestrial reference systems, but can also be used to determine the very-long wavelength part of the Earth's gravity field.

Leakage: Leakage describes the inability to localize signals in the GRACE-derived data sets exactly, due to, among other factors, the band-limited resolution of GRACE and the filtering. Together, leakage leads to apparent signal loss (leakage out) or gain (leakage in) inside a given integration region.

Level-0 data: Raw instrument data received from the satellites. 

Level-1 data: Level-0 data converted to engineering units (Level-1A). Time-tagging of the data of all different instruments, filtering, and downsampling results in Level-1B data. Level-1B data includes range rate data, GPS positions, and accelerometer data.

Level-2 data: Time-variable, usually monthly, gravity fields given in spherical harmonic (SH) coefficients.

Level-2B data: Post-processed Level-2 data, still given in SH coefficients. The post-processing includes several corrections that are applied to obtain as accurate as possible TWSA data.

Level-3 data: Time-variable, usually monthly, gravity fields given in the spatial domain as gridded data, based on Level-2B data. TWSA data is a Level-3 data set. 

Executive Summary

This document serves as the Algorithm Theoretical Basis Document (ATBD) for the Copernicus Climate Change Service (C3S) Terrestrial Water Storage Anomaly (TWSA) data product (version v1.0), derived from observations made by the Gravity Recovery and Climate Experiment (GRACE) and GRACE-Follow-On (FO) satellite missions.

The ATBD outlines the processing steps required to generate gridded Level-3 TWSA data from Level-2 GRACE and GRACE-FO data. The TWSA production system picks up Level-2 data from the Combination Service for Time-variable Gravity Fields (COST-G) service. Accordingly, this document describes the processing steps beyond Level-2. Additionally, it details all necessary datasets, including input, external, and auxiliary data.

Section 1 introduces the GRACE and GRACE-FO satellite missions and explains their measurement principles in some detail.

Section 2 presents all datasets used in the production system. Section 2.1 describes the Level-2 GRACE and GRACE-FO data processed externally, outside the C3S production system. Section 2.2 summarises the auxiliary datasets required in the production system, including low-degree spherical harmonic (SH) coefficients, glacial isostatic adjustment (GIA) correction, and the land-ocean mask.

Section 3 details the production system, dividing the processing into two stages: spherical harmonic (SH) domain processing (Level-2B data, Section 3.1) and spatial domain processing (Level-3 data, Section 3.2). Level-2B processing includes the replacement of specific low-degree SH coefficients and GIA correction. Level-3 processing involves transforming data from the SH domain to the spatial (gridded) domain and correcting for co- and post-seismic signals from megathrust earthquakes. Section 3.2 also explains the uncertainty estimation of the final TWSA product.

Finally, Section 4 summarises the output fields of the final TWSA product.

Missions and Instruments

Table 1: Details of missions. Operation dates refer to the periods where full sets of parameters for gravity field calculation could be obtained.

MissionStart date of operationEnd date of operationOrbit altitudeOrbit inclinationPrimary instruments
GRACE04/04/200229/06/2017
  • Launch: 500 km
  • End of operation:  330 km

89°

  • K/Ka-band ranging assembly
  • GPS receiver
  • ONERA SuperSTAR accelerometer
  • Star camera assembly
GRACE-FO01/06/2018in operation
  • Launch: 490 km
89°
  • K/Ka-band ranging assembly
  • GPS receiver
  • ONERA SuperSTAR accelerometer
  • Star camera assembly
  • Laser Ranging Interferometer (serving as technology demonstration)

Terrestrial Water Storage Anomaly (TWSA) is derived solely from data obsevered by the Gravity Recovery and Climate Experiment (GRACE) (2002–2017, Tapley et al., 2004) and GRACE-Follow-On (FO) (since 2018, Landerer et al., 2020) satellite missions (details see Table 1). GRACE was a joint US-German mission led by the National Aeronautics and Space Administration (NASA) and the German Aerospace Center (DLR). GRACE-FO is also a US-German collaboration, this time between NASA and the Helmholtz Centre for Geosciences (GFZ).

Both missions consist of two satellites in the same polar orbit, separated by approximately 220 km. A microwave interferometer — supplemented by an additional laser interferometer on GRACE-FO — continuously measures the distance between the satellites and its variations. The absolute position and orientation of the satellites are determined using Global Positioning System (GPS) receivers and star cameras, while accelerometers measure non-gravitational forces acting on the satellites. Unequal density distributions within the Earth's interior, along with spatially varying surface masses, lead to a spatially variable gravity field. These mass variations influence the orbital velocities of the satellites.

Figure 1 illustrates the measurement principle of GRACE and GRACE-FO. Both satellites experience the same gravity field anomalies along their orbit, but with a slight time delay of approximately 30 seconds. As a result, their relative velocity changes continuously, causing small variations in the distance between them, which are constantly measured. Typically, data are collected over a month and, in combination with measurements from other onboard instruments, used to derive a monthly gravity field model. The GRACE and GRACE-FO missions have thus enabled, for the first time, the observation of monthly, seasonal, and longer-term mass variations in the Earth's system. These variations are primarily attributed to changes in surface mass distribution, particularly water mass redistribution.


Figure 1: Infographic illustrating the measuring method of the GRACE and GRACE-FO satellites. (1) The distance between the two satellites remains nearly constant over a surface with little or no spatial mass variation. (2) As the first satellite approaches a mass anomaly (e.g., a mountain), it experiences an acceleration due to the increased gravitational pull. However, the second satellite, following 220 km behind, is not yet affected and continues at its original velocity. As a result, the distance between the satellites increases. (3) Once the first satellite moves past the mass anomaly, it slows down. Meanwhile, the second satellite, now reaching the anomaly, experiences acceleration in the same way the first satellite did earlier. This causes another change in the distance between them. (4) Finally, after the trailing satellite has also passed the mass anomaly and slowed to its original velocity, the initial state is restored, and the distance between the two satellites returns to nearly constant again. © Klinghammer, GFZ


Input and Auxiliary Data

GRACE and GRACE-FO Data

GRACE and GRACE-FO data processing is categorised into different science data levels. First, raw instrument data (Level-0) received from the satellites are converted into engineering units (Level-1A). Next, Level-1A data are processed into Level-1B data through steps such as uniform time-tagging of individual data products from different instruments, filtering, and down-sampling of raw data rates. From Level-1B data, time-variable gravity field models (Level-2) are derived.

This gravity field retrieval is conducted by several institutes worldwide, including the three official Science Data System (SDS) processing centres for GRACE and GRACE-FO Level-2 data: GFZ, the Center for Space Research at the University of Texas (CSR), and NASA’s Jet Propulsion Laboratory (JPL). Additional unofficial Level-2 solutions are provided by several institutions, e.g. the Astronomical Institute of the University of Bern (AIUB), Switzerland, the Institute of Geodesy at Graz University of Technology (TUG), Austria, or HuaZhong University of Science and Technology (HUST), People's Republic of China. While all these institutes use similar methods, each employs a slightly different approach. Further details on current Level-2 processing can be found, for example, in Dahle et al. (2019a) for GFZ, Meyer et al. (2016) for AIUB, or Kvas et al. (2019) for TUG. 

To provide users with the best monthly global gravity field model available, the Combination Service for Time-variable Gravity Fields (COST-G, URL last accessed 17/03/2025) was established (Jäggi et al., 2020). COST-G operates as a product centre of the International Gravity Field Service (IGFS) under the umbrella of the International Association of Geodesy (IAG). It combines monthly global gravity field models from the above-mentioned Level-2 processing centres (Meyer et al., 2024). These combined gravity field models benefit from the strengths of the various processing approaches of the contributing centres, resulting in a reduced noise level. The C3S TWSA product is based on the latest GRACE and GRACE-FO data products provided by COST-G RL02. Typically, new COST-G monthly gravity models are made publicly available with a latency of two to three months.

Due to mission operations or instrument issues, some monthly solutions are not available or have lower data quality. The tables 2.1.1., 2.1.2, and the Product User Guide and Specification (PUGS) detail the available months.

COST-G (GRACE)

Originating System

Gravity Recovery and Climate Experiment (GRACE)

Data Class

Earth observation

Key Technical Characteristics

  • Inter-satellite K/Ka-Band ranging (KBR) system, precision 10 μm
  • Orbit height: 500 km (at launch, declining orbit); inclination angle: 89°; period of revolution: ~90 min; inter-satellite distance: ~220 km
  • Provided as SH coefficients with maximum degree/order 90, related to the spatial wavelength of 220 km, effective spatial resolution ~300 km
  • Temporal resolution: monthly

Data Availability and Coverage

April 2002 - June 2017, global

Source Data Name and Product Technical Specifications

COST-G RL02

Data set reference:

  • Meyer, U, Jäggi, A, Dahle, C, Flechtner, F, Kvas, A, Behzadpour, S, Öhlinger, F, Mayer-Gürr, T, Lemoine, J-M, Bourgogne, S, Lasser, M, Koch, I, Flury, J, Chen, Q, Wang, C, Yan, Z, Zhou, H, and Feng, W 2025 International Combination Service for Time-variable Gravity Fields (COST-G) Monthly GRACE/GRACE-FO RL02 Series. DOI: https://doi.org/10.5880/COST-G.ICGEM_02_L2

Technical Specification:

  • Meyer et al. (2024) 

Data set:

Data Quantity

Total volume is 22.4 MB (compressed)

Data Quality and Reliability

Instrument specification:

  • Tapley et al. (2004)

Validation report:

  • Jäggi et al. (2020)
  • Meyer et al. (2024)

Ordering and Delivery Mechanism

No ordering necessary

Delivery via the International Centre for Global Earth Models (ICGEM):

Access Conditions and Pricing

Freely accessible

Issues

  • Missing monthly solutions due to unavailability of all required satellite data:
    • 06/2002, 07/2002, 06/2003, 01/2011, 06/2011, 05/2012, 10/2012, 03/2013, 08/2013, 09/2013, 02/2014, 07/2014, 12/2014, 06/2015, 10/2015, 11/2015, 04/2016, 09/2016, 10/2016, 02/2017
  • Months affected by short repeat orbits leading to solutions with poorer data quality, higher uncertainty, and larger data latency due to more complex data processing:
    • 09/2004, affected solutions: 08/2004, 09/2004, 10/2004
    • 05/2012, affected solutions: 04/2012, 06/2012
    • 12/2013, affected solutions: 12/2013, 01/2014
    • 02/2015, affected solutions: 01/2015, 02/2015
  • Starting with the solution for 11/2016, no accelerometer data on GRACE-B is available. The data has to be transplanted from the accelerometer data of GRACE-A (details: Bandikova et al., 2019).

COST-G (GRACE-FO)

Originating System

GRACE Follow-On (GRACE-FO)

Data Class

Earth observation

Key Technical Characteristics

  • Inter-satellite K/Ka-Band ranging (KBR) system, precision 10 μm; Laser Ranging Interferometer, precision 10nm
  • Orbit height: 490 km (at launch, declining orbit); inclination angle: 89°; period of revolution: ~90 min; inter-satellite distance: ~220 km
  • Provided as SH coefficients with maximum degree/order 90, related to the spatial wavelength of 220 km, effective spatial resolution ~300 km
  • Temporal resolution: monthly

Data Availability and Coverage

June 2018 - on going, global

Source Data Name and Product Technical Specifications

COST-G RL02

Data set reference:

  • Meyer, U, Jäggi, A, Dahle, C, Flechtner, F, Kvas, A, Behzadpour, S, Öhlinger, F, Mayer-Gürr, T, Lemoine, J-M, Bourgogne, S, Lasser, M, Koch, I, Flury, J, Chen, Q, Wang, C, Yan, Z, Zhou, H, and Feng, W 2025 International Combination Service for Time-variable Gravity Fields (COST-G) Monthly GRACE/GRACE-FO RL02 Series. DOI: https://doi.org/10.5880/COST-G.ICGEM_02_L2

Technical Specification:

  • Meyer et al. (2024)

Data set:

Data Quantity

Total volume as of 18/03/2025: 10.6 MB (compressed )

Data Quality and Reliability

Instrument specification:

  • Landerer et al. (2020)

Validation report:

  • Jäggi et al. (2020)
  • Meyer et al. (2024)

Ordering and Delivery Mechanism

No ordering necessary

Delivery via the International Centre for Global Earth Models (ICGEM)

Access Conditions and Pricing

Freely accessible

Issues

  • Missing monthly solutions due to unavailability of all required satellite data:
    • 08/2018, 09/2018
  • Months affected by short repeat orbits leading to solutions with poorer data quality, higher uncertainty, and larger data latency due to more complex data processing:
    • 04/2023, affected solutions: 04/2023, 05/2023

    • 09/2024, affected solutions: 09/2024, 10/2024

  • Accelerometer data on the second GRACE-FO satellite (GRACE-D) is below expected quality since launch. While GRACE-D accelerometer data can partially be used, transplanted data from the accelerometer of the first GRACE-FO satellite (GRACE-C) is additionally needed (details: Harvey et al., 2024).

Auxiliary Data

Low-Degree Spherical Harmonic Coefficients

It is recommended that specific spherical harmonic (SH) coefficients estimated from GRACE and GRACE-FO are replaced with more reliable external estimates. In particular, the SH coefficient of degree 2 and order 0 (C20), which relates to the Earth's flattening, is affected by systematic errors (Cheng and Ries, 2017). Depending on the quality of accelerometer observations, the C30 coefficient is poorly determined, too, during specific periods (Loomis et al., 2020). Commonly, C20 and C30 coefficients derived from satellite laser ranging (SLR) observations of dedicated SLR satellites or a combination of GRACE, GRACE-FO, and SLR data are used for this recommended replacement.

For the TWSA production system, the latter (GRACE, GRACE-FO, and SLR) data are used and processed at GFZ. The SLR processing follows the same methodology as GFZ’s SLR-only C20 time series (König et al., 2019), incorporating observations from six satellites: Laser Geodynamics Satellite 1 (LAGEOS-1), LAGEOS-2, Ajisai, Stella, Satellite de Taille Adaptée avec Réflecteurs Laser pour les Études de la Terre (Starlette), and Laser Relativity Satellite (LARES). The combination with GRACE and GRACE-FO is performed at the level of normal equations, using relative weights assigned to individual SLR satellites based on variance component estimation. Additionally, an empirically derived relative weighting is applied between SLR and GRACE and GRACE-FO. From the combined normal equations, gravity field solutions up to degree and order 6 are estimated independently for each month, from which the SH coefficients C20 and C30 are extracted.

Originating System

Passive SLR satellites (LAGEOS-1, LAGEOS-2, Ajisai, Stella, Starlette, LARES), GRACE, GRACE-FO

Data Class

Earth observation

Key Technical Characteristics

For SLR:

For GRACE and GRACE-FO:

  • see tables in sections 2.1.1 and 2.1.2

For the combination of SLR and GRACE/GRACE-FO:

  • Spatial resolution: maximum SH degree and order 6
  • Temporal resolution: monthly

Data Availability and Coverage

April 2002 - on going, global

Source Data Name and Product Technical Specifications

Technical Specification:

  • SLR: König et al. (2019)
  • GRACE/GRACE-FO: Dahle et al. (2019a)

Data sets:

  • SLR: König et al. (2019)
  • GRACE: Dahle et al. (2018)
  • GRACE-FO: Dahle et al. (2019b)

Data Quantity

Total volume as of 18/03/2025: <1 MB 

Data Quality and Reliability

Data is already successfully used for GFZ's operationally generated GRACE/GRACE-FO Level-3 products (Dahle et al., 2025) 

Ordering and Delivery Mechanism

not applicable (in-house processing at GFZ)

Access Conditions and Pricing

not applicable (in-house processing at GFZ)

Issues

  • SLR processing depends on the availability of data products provided by the International Laser Ranging Service (ILRS)

Model of Glacial Isostatic Adjustment

Glacial Isostatic Adjustment (GIA) refers to the surface deformation caused by ice-mass changes over the last 100,000 years, dominated by the end of the most recent glacial cycle. Due to the Earth's viscoelastic response to the mass redistribution between the grounded ice sheets and the fluid ocean, the Earth's gravity field is affected by secular trends, particularly in formerly glaciated regions in North America, Fennoscandia, and Antarctica. 

To account for GIA effects, the ICE-6G_D (VM5a) model (Peltier et al., 2018) is used to correct the spherical harmonic coefficients. Figure 2 below presents the GIA uplift rates converted to a global 1° x 1° grid.

Figure 2: Annual GIA uplift rates from the applied ICE-6G_D correction, expressed in equivalent water height change per year.


Data class

Model

Key technical characteristics

  • GIA model developed to fit GPS observations from North America, Eurasia, and Antarctica
  • Improved bathymetry for the Southern Ocean

Data Availability and Coverage

Given in SH coefficients up to degree and order 256

Source Data Name and Product Technical Specifications

ICE-6G_D (VM5a)

Technical Specification: Peltier et al. (2018)

Data Quantity

1.5 MB

Data Quality and Reliability

not applicable

Ordering and delivery mechanism

Download from the Supporting Information of Peltier et al. (2018)

https://agupubs.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1002%2F2016JB013844&file=jgrb52450-sup-0003-Data_S3.txt (URL last accessed 18/03/2025)

Access conditions and pricing

Freely accessible 

Issues

No issues known

Land-Ocean Mask

For the TWSA data set, a global land-ocean mask is required to distinguish continental areas from the oceans. The mask used in the TWSA production system is part of a data set produced during the G3P project based on the European Space Agency (ESA) Climate Change Initiative (CCI) global map of open water bodies v4.0 (CCI WB v4.0) and are freely available in different resolutions (Sharifi et al., 2025). Here we consider all grid cells of the TWSA data set (0.5° x 0.5° resolution) as land if they contain any land. This is illustrated for Europe in Figure 3.

Figure 3: Land-ocean mask used in the TWSA product for Europe (land: green, ocean: blue), overlaid by the coastline of the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG, Wessel and Smith, 1996). 

Originating System

Multiple individual synthetic aperture radar (SAR) and optical water body and auxiliary datasets, for full list see Lamarche et al. (2017), section 2, table 2.

Data class

Earth observation

Key technical characteristics

  • Static global map of open water bodies

  • Version used in this product distinguishes between ocean and inland. All inland water bodies are considered land.

  • Spatial resolution: 150 m

  • Legend: 0-Ocean, 1-Land

  • Available as tif

Data Availability and Coverage

Void-free full global coverage

Source Data Name and Product Technical Specifications

European Space Agency (ESA) Climate Change Initiative (CCI) global map of open water bodies v4.0 (CCI WB v4.0)
Technical Specification:
Lamarche et al. (2017)

Data Quantity

298 MB (original data), 4 MB (G3P land-ocean masks)

Data Quality and Reliability

Validated against 2110 samples. Overall accuracy between 98% and 100%. Class representation F-score of 89%.
Validation report in
Lamarche et al., 2017, section 3.

Ordering and delivery mechanism

G3P land-ocean masks available from the GFZ data service:
https://doi.org/10.5880/G3P.2025.001

Original data set downloadable after registration from ESA CCI data viewer:
https://maps.elie.ucl.ac.be/CCI/viewer/download.php
(URL last accessed 18/03/2025)

Access conditions and pricing

Freely accessible

Issues

No issues known

The CCI WB mask is resampled to a global 0.5° x 0.5° grid. Hereby, each grid cell containing a land cell in the original CCI WB 150m resolution is categorised as land. The resulting shoreline is illustrated for Europe in Figure 3.

Algorithms

This section provides information about the algorithms of the TWSA production system. The reader is referred to Dahle et al. (2025) for more details on the algorithms. The key differences between the TWSA product presented here and the product described in Dahle et al. (2025) are highlighted in Sections 3.1 and 3.2.

Figure 4 provides an overview of the data processing workflow, from Level-2 (COST-G spherical harmonic coefficients) to Level-3 (gridded TWSA data).


Figure 4: Overview of the processing steps from GRACE and GRACE-FO Level-2 data to the TWSA (Level-3) data product. The processing steps leading to COST-G Level-2 data are explained in section 2.1 and not part of the C3S TWSA production system. The processing from Level-2 to Level-2B is conducted in the spherical harmonic domain (indicated in green), whereas the subsequent processing steps from Level-2B to the final TWSA product are performed in the spatial domain (indicated in blue). Additional independent input data required for the processing of TWSA are referred to as auxiliary data.

Level-2B

Reduction of Long-Term Mean

A long-term mean is subtracted from the monthly gravity fields to obtain anomalies. Here, the mean gravity field is computed using all available GRACE and GRACE-FO monthly fields from January 2003 to December 2022 (see PUGS and Sections 2.1.1 and 2.1.2 for details on the available months). 

Anisotropic Filtering

The observation geometry primarily depends on inter-satellite ranging in the along-track direction at near-polar orbits. This results in spatially correlated errors in the GRACE and GRACE-FO gravity field solutions, leading to anisotropic error patterns in the north-south direction, known as striping errors. To mitigate these errors and effectively separate signal from noise in the GRACE and GRACE-FO Level-2 products, filtering is required.

The widely used decorrelation method by Kusche et al. (2009), known as the DDK filter, was adapted by Horvath et al. (2018) to consider temporal variations in filter strength based on error covariances. This adapted filter, the time-variable DDK (VDK) filter, is used in the TWSA production system. The strength of VDK filters is indicated by a trailing number, where 1 represents the strongest commonly available filter. In the TWSA production system, the SH coefficients are filtered using VDK2, VDK3, and VDK5. Figure 5 illustrates the effect of these three different filter strengths on the March 2008 gravity field. More details on anisotropic filtering and the advantages of VDK filtering above DDK filtering can be found in Section 2.1.1 and Figure 2 in Dahle et al. (2025).

All subsequent processing steps are carried out in parallel for the three sets of coefficients.


Figure 5: Comparison between the three different filter strengths for the same month, March 2008: (a) shows the VDK2 filter which has strongest smoothing of the filters used in the TWSA production system; (b) shows the intermediate smoothing VDK3 filter where slight stripping residuals are visible in the Pacific ocean; (c) shows the weakest VDK5 filter of the TWSA production system where a significant amount of striping residuals remained.

The VDK5 filter retains the signal amplitude and localisation more clearly. On the other hand, the VDK2 filter smoothes the signals.

Replacement of Specific Low-Degree Spherical Harmonic Coefficients

As explained in Section 2.2.1, specific low-degree SH coefficients (C20 for the entire time series and C30 from November 2016 to June 2017) are replaced by estimates derived from a combination of SLR, GRACE, and GRACE-FO observations. More details about replacing the low-degree SH coefficients and the consequences of the chosen replacement data set can be found in Section 2.1.2 in Dahle et al. (2025).

Geophysical Correction of Signals Induced by Glacial Isostatic Adjustment

The GIA model ICE-6G-D (see Section 2.2.2) is subtracted from the SH coefficients to correct for the mass effect of GIA uplift. Section 2.1.5 in Dahle et al. (2025) provides more information on the GIA correction.

Approximation of Geocentre Variations

The SH coefficients of degree 1 (C10, C11, S11) are related to the distance between the Earth’s centre of figure (CF) and centre of mass (CM), commonly referred to as geocentre motion. Since the GRACE and GRACE-FO satellite missions are insensitive to CF, these degree-1 SH coefficients are not estimated and are instead set to zero by definition.

Nevertheless, this information is essential for accurately quantifying oceanic and terrestrial mass distributions. To incorporate it into the Level-2B products, external estimates of degree-1 coefficients are required. The TWSA production system applies the approximation method by Swenson et al. (2008) to derive degree-1 coefficients. This approximation is based on monthly GRACE and GRACE-FO SH coefficients, and the estimated degree-1 values are then inserted into the monthly Level-2B products. Section 2.1.3 in Dahle et al. (2025) details the approximation scheme used here.

Empirical Correction of S2 Tidal Aliasing Errors

Removing ocean tidal signals from satellite observations is essential during GRACE and GRACE-FO Level-2 processing. However, temporal aliasing of unmodelled or insufficiently modelled ocean tide signals caused by the space-time sampling of a satellite gravimetry mission, introduces additional gravity field errors. A particularly significant aliasing frequency in the GRACE and GRACE-FO gravity fields has a 161-day period, which may be linked to model errors in the semi-diurnal solar tide S2 present in both the ocean and atmosphere. More background information on this S2 tidal aliasing error can be found in Section 2.1.4. in Dahle et al. (2025).

To mitigate these S2 tidal aliasing errors, an empirical correction is implemented into in the TWSA production system. This involves simultaneously fitting a bias, linear trend, annual, semi-annual, and 161-day periodic signals to the time series of monthly spherical harmonic (SH) coefficients. The 161-day periodic signal is then evaluated at the mid-epoch of each month and subsequently subtracted.

The deterministic signal components are estimated only once, based on the data used to compute the long-term mean (see Section 3.1.1), and are subsequently extrapolated. A phase offset of 100° is applied between GRACE and GRACE-FO when fitting the S2 tidal aliasing frequency, as the nodal planes of their orbits are not aligned (Landerer et al., 2020).

Level-3

Spherical Harmonic Synthesis

The spherical harmonic synthesis transforms the Level-2B SH coefficients to globally gridded data on a 0.5°x0.5° longitude-latitude grid (Wahr et al., 1998). The grid used in this product differs from that described in Dahle et al. (2025). The spherical harmonic synthesis is described by equation (1):

\[ twsa(\lambda, \theta) [mm] = 1000\frac{1}{\rho_W} \frac{M}{2 \pi R^2} \sum_{n=0}^{N} \frac{2n+1}{1+k'_n} \sum_{m=0}^{n}(c_{nm}P_{nm}(\cos \theta) \cos m\lambda + s_{nm}P_{nm}(\cos \theta) \sin m\lambda) \quad (1)\]

The gridded data is expressed in equivalent water heights in millimetres at a given location with longitude λ and co-latitude θ. cnm and snm are the spherical harmonic coefficients of degree n and order m; Pnm the fully normalised Legendre functions; kn'  the Load-Love-Numbers; ρw density of water; and M and R mass and radius of the Earth. The spherical harmonic synthesis is done up to the maximum degree and order N. The reference surface for the spherical harmonic synthesis to the grid is the reference ellipsoid defined in the International Earth Rotation and Reference Systems Service (IERS) Conventions (2010) Tab 1.1 (Petit and Luzum, 2010). 

The spherical harmonic synthesis is performed for Level-2B data filtered with VDK2, VDK3, and VDK5.

Combination of Different Filter Strengths

All gridded datasets (filtered using VDK2, VDK3, and VDK5) are decomposed into a linear trend component, annual and semi-annual harmonics, and residuals. The residuals contain not only signal noise, including residual striping errors that persist after filtering, but also inter-annual variability. Hence, deterministic signals can be reliably estimated from the weaker VDK5 filtered fields while being less affected by leakage, providing higher spatial resolution. Therefore, in the TWSA production system, the deterministic signals from VDK5-filtered fields are combined with the residual signals from VDK3-filtered fields to minimise leakage effects in the final dataset.

Prior to this combination, the noise level of the residual signal from the VDK3-filtered field is evaluated against the expected noise level of the entire time series. If, in a given month, the standard deviation of residuals in the open ocean (distance to coast > 1000 km) exceeds twice the mean of all monthly standard deviations, the residual variations in the VDK3-filtered grids are replaced with those from the VDK2-filtered grids. The TWSA dataset includes a flag indicating the filter strength used in each case.

Geophysical Correction of Co- and Post-Seismic Deformations from Megathrust Earthquakes

Megathrust earthquakes generate gravity signals by displacing large volumes of tectonic plate material. These mass changes must be removed from the TWSA product for hydrological applications. The three earthquakes considered in the TWSA production system were  (1) Sumatra–Andaman in 2004, (2) Maule, Chile, in 2010, and (3) Tohoku-oki, Japan, in 2011. For each event, the position and timing are obtained from an earthquake catalogue. A step function is then fitted within a spherical cap (1000 km radius) centred at the epicentre to all gridded monthly solutions. An exponential decay function is applied over two years following the main event. These empirical estimates of co- and post-seismic gravity signals are then subtracted from the monthly gravity field solutions, ensuring that only hydrology-related signals are retained. More information on the background of the earthquake correction can be found in Section 2.1.6 in Dahle et al. (2025).

Uncertainty Estimation

The precision and accuracy of the TWSA data grids are jointly described by uncertainty expressed as a standard deviation. These uncertainties are influenced, among other factors, by the uncertainties of the sensors on the satellites, uncertainties of background models applied in the processing chain, the orbit configuration of the missions (near-polar orbit), or environmental effects such as solar activity.

The uncertainty estimation of TWSA is based on the open ocean (distance to coast >1000km) area weighted standard deviation σ(t) for each month t defined in equation (2), after removing the deterministic signals, i.e. trend, annual, and semiannual signal (twsared). 

\[\sigma(t) = \sqrt{\frac{1}{\sum_i^a w_i}\sum_{i}^a w_i (twsa_{red}(\lambda_i, \theta_i,t)-\overline{twsa_{red}(t)})^2} \quad (2)\]

With a representing the number of open ocean grid points, and wi being the area weight of the ith grid cell at the longitude λi and latitude θi. This monthly value is latitude-dependently scaled according to the covariance model theoretically described in Section 3 in Boergens et al. (2020). Table 1 in Boergens et al. (2022) provides the numerical values of the model parameters. The covariance model describes the covariance between two points as given in equation (1) in Boergens et al. (2020). For the latitude-depending scaling of σ(t) the corner case of two identical points is assumed.

This uncertainty estimation is not included in the data product of Dahle et al. (2025).

Land-Ocean Masking

The final step of the TWSA data product processing is masking out the oceans using the land-ocean mask described in Section 2.2.3. Please note that the land-ocean mask applied in the TWSA product differs from the one used by Dahle et al. (2025).

Output Data

Table 2 lists the output fields generated by the TWSA production system and made available to users. Detailed information on each field and the file format are provided in the Product User Guide and Specification (PUGS).

Table 2: Overview of output data fields.

Field

Description

Unit

twsa

Gravity-based terrestrial water storage anomaly

mm

twsa_uncertainty

Uncertainty of the terrestrial water storage anomaly

mm

flag_filter

Flag indicating the filter strength (VDKx) used for the residual interannual signal in each month.

For example: 2=VDK2 used, 3=VDK3 used

-

Figure 6 shows the final data product for March 2008 with (a) TWSA and (b) TWSA uncertainty. In (c) the flag_filter field is shown for the entire time series. Currently, only for two months - February 2015 and December 2016 - the residuals require filtering with the stronger VDK2 filter.

 

Figure 6: Fields of the final data product. (a) shows TWSA from the final data product for the exemple month of March 2008 stored in the twsa variable in the netcdf file; (b) shows the uncertainty of TWSA for the same month, stored in the twsa_uncertainty variable of the netcdf file; (c) displays the time series of the used filter strength for the interannual signal which is stored in the flag_filter variable of the netcdf file.

References

  • Bandikova, T., McCullough, C., Kruizinga, G.L., Save, H., Christophe, B., 2019. GRACE accelerometer data transplant. Advances in Space Research 64, 623–644. https://doi.org/10.1016/j.asr.2019.05.021  [URL last accessed 17/03/2025]

  • Boergens, E., Dobslaw, H., Dill, R., Thomas, M., Dahle, C., Murböck, M., Flechtner, F., 2020. Modelling spatial covariances for terrestrial water storage variations verified with synthetic GRACE-FO data. Int J Geomath 11, 24. https://doi.org/10.1007/s13137-020-00160-0  [URL last accessed 17/03/2025]

  • Boergens, E., Kvas, A., Eicker, A., Dobslaw, H., Schawohl, L., Dahle, C., Murböck, M., Flechtner, F., 2022. Uncertainties of GRACE-Based Terrestrial Water Storage Anomalies for Arbitrary Averaging Regions. Journal of Geophysical Research: Oceans. https://doi.org/10.1029/2021JB022081 [URL last accessed 17/03/2025]
  • Cheng, M., Ries, J., 2017. The unexpected signal in GRACE estimates of C20. J. Geodesy, 91, 897–914. https://doi.org/10.1007/s00190-016-0995-5 [URL last accessed 17/03/2025]
  • Dahle, C., Flechtner, F., Murböck, M., Michalak, G., Neumayer, H., Abrykosov, O., Reinhold, A., König, R., 2018. GRACE Geopotential GSM Coefficients GFZ RL06, V. 6.0. GFZ Data Services. https://doi.org/10.5880/GFZ.GRACE_06_GSM [URL last accessed 17/03/2025]

  • Dahle, C., Murböck, M., Flechtner, F., Dobslaw, H., Michalak, G., Neumayer, K.H., Abrykosov, O., Reinhold, A., König, R., Sulzbach, R., Förste, C., 2019a. The GFZ GRACE RL06 monthly gravity field time series: Processing details and quality assessment. Remote Sensing 11(18), 2116. https://doi.org/10.3390/rs11182116 [URL last accessed 17/03/2025]
  • Dahle, C., Flechtner, F., Murböck, M., Michalak, G., Neumayer, K. H., Abrykosov, O., Reinhold, A., König, R., 2019b. GRACE-FO Geopotential GSM Coefficients GFZ RL06, V. 6.3. GFZ Data Services. https://doi.org/10.5880/GFZ.GRACEFO_06_GSM [URL last accessed 17/03/2025]

  • Dahle, C., Boergens, E., Sasgen, I., Döhne, T., Reißland, S., Dobslaw, H., Klemann, V., Murböck, M., König, R., Dill, R., Sips, M., Sylla, U., Groh, A., Horwath, M., Flechtner, F., 2025. GravIS: mass anomaly products from satellite gravimetry. Earth System Science Data 17, 611–631. https://doi.org/10.5194/essd-17-611-2025 [URL last accessed 17/03/2025]
  • Harvey, N., Bertiger, W., McCullough, C., Miller, M., Save, H., Yuan, D.-N., 2024. Recovering differential forces from the GRACE-D accelerometer. Earth Space Sci., 11, e2023EA003200. https://doi.org/10.1029/2023EA003200 [URL last accessed 17/03/2025]
  • Horvath, A., Murböck, M., Pail, R., Horwath, M., 2018. Decorrelation of GRACE time variable gravity field solutions using full covariance information. Geosciences 8, 323. https://doi.org/10.3390/geosciences8090323 [URL last accessed 17/03/2025]
  • Jäggi, A., Meyer, U., Lasser, M., Jenny, B., Lopez, T., Flechtner, F., Dahle, C., Förste, C., Mayer-Gürr, T., Kvas, A., Lemoine, J.-M., Bourgogne, S., Weigelt, M., Groh, A., 2020. International Combination Service for Time-Variable Gravity Fields (COST-G): Start of Operational Phase and Future Perspectives, in: International Association of Geodesy Symposia. Springer Berlin Heidelberg, Berlin, Heidelberg. https://doi.org/10.1007/1345_2020_109 [URL last accessed 17/03/2025]
  • König, R., Schreiner, P., Dahle, C., 2019. Monthly estimates of C(2,0) generated by GFZ from SLR satellites based on GFZ GRACE and GRACE-FO RL06 background models, V. 1.0. GFZ Data Services. https://doi.org/10.5880/GFZ.GRAVIS_06_C20_SLR [URL last accessed 17/03/2025]
  • Kusche, J., Schmidt, R., Petrovic, S., Rietbroek, R., 2009. Decorrelated GRACE time-variable gravity solutions by GFZ, and their validation using a hydrological model. J Geod 83, 903–913. https://doi.org/10.1007/s00190-009-0308-3 [URL last accessed 17/03/2025]

  • Kvas, A., Behzadpour, S., Ellmer, M., Klinger, B., Strasser, S., Zehentner, N., Mayer-Gürr, T., 2019. ITSG-Grace2018: Overview and Evaluation of a New GRACE-Only Gravity Field Time Series. J. Geophys. Res.-Sol. Ea. 124, 9332–9344. https://doi.org/10.1029/2019JB017415 [URL last accessed 17/03/2025]
  • Landerer, F.W., Flechtner, F.M., Save, H., Webb, F.H., Bandikova, T., Bertiger, W.I., Bettadpur, S.V., Byun, S., Dahle, C., Dobslaw, H., Fahnestock, E., Harvey, N., Kang, Z., Kruizinga, G.L.H., Loomis, B.D., McCullough, C., Murböck, M., Nagel, P., Paik, M., Pie, N., Poole, S., Strekalov, D., Tamisiea, M.E., Wang, F., Watkins, M.M., Wen, H. ‐Y., Wiese, D.N., Yuan, D. ‐N., 2020. Extending the global mass change data record: GRACE Follow‐On instrument and science data performance. Geophys. Res. Lett. https://doi.org/10.1029/2020GL088306 [URL last accessed 17/03/2025]
  • Lamarche, C., Santoro, M., Bontemps, S., D’Andrimont, R., Radoux, J., Giustarini, L., Brockmann, C., Wevers, J., Defourny, P., Arino, O., 2017. Compilation and Validation of SAR and Optical Data Products for a Complete and Global Map of Inland/Ocean Water Tailored to the Climate Modeling Community. Remote Sensing 9, 36. https://doi.org/10.3390/rs9010036 [URL last accessed 17/03/2025]

  • Loomis, B. D., Rachlin, K. E., Wiese, D. N., Landerer, F. W., Luthcke, S. B., 2020. Replacing GRACE and GRACE-FO C30 with satellite laser ranging: Impacts on Antarctic Ice Sheet mass change. Geophys. Res. Lett., 47, e2019GL085488. https://doi.org/10.1029/2019GL085488 [URL last accessed 17/03/2025]
  • Meyer, U., Jäggi, A., Jean, Y., Beutler, G., 2016. AIUB-RL02: an improved time-series of monthly gravity fields from GRACE data. Geophys. J. Int. 205, 1196–1207. https://doi.org/10.1093/gji/ggw081 [URL last accessed 17/03/2025]
  • Meyer, U., Lasser, M., Dahle, C., Förste, C., Behzadpour, S., Koch, I., Jäggi, A., 2024. Combined monthly GRACE-FO gravity fields for a Global Gravity-based Groundwater Product. Geophys. J. Int., 236, 456–469. https:/doi.org/10.1093/gji/ggad437 [URL last accessed 17/03/2025]
  • Peltier, R.W., Argus, D.F., Drummond, R., 2018. Comment on “An Assessment of the ICE-6G_C (VM5a) Glacial Isostatic Adjustment Model” by Purcell et al.: The ICE-6G_C (VM5a) GIA model. J. Geophys. Res. Solid Earth 123, 2019–2028. https://doi.org/10.1002/2016JB013844 [URL last accessed 17/03/2025]
  • Petit, G., Luzum, B., 2010. IERS Conventions 2010 (No. 36), IERS Technical Note. International Earth Rotation and Reference Systems Service, Frankfurt.
  • Sharifi, E., Haas, J., Boergens, E., Dobslaw, H., and Güntner, A., 2025a. Technical note: GRACE-compatible filtering of water storage data sets via spatial autocorrelation analysis, Hydrol. Earth Syst. Sci., 29, 6985–6998, https://doi.org/10.5194/hess-29-6985-2025.
  • Tapley, B.D., Bettadpur, S., Watkins, M., Reigber, C., 2004. The gravity recovery and climate experiment: Mission overview and early results. Geophys. Res. Lett. 31. https://doi.org/10.1029/2004GL019920 [URL last accessed 17/03/2025]
  • Wahr, J., Molenaar, M., Bryan, F., 1998. Time variability of the Earth’s gravity field: Hydrological and oceanic effects and their possible detection using GRACE. Journal of Geophysical Research 103, 30205–30230.
  • Wessel, P., Smith, W.H.F., 1996. A global, self-consistent, hierarchical, high-resolution shoreline database. Journal of Geophysical Research: Solid Earth 101, 8741–8743. https://doi.org/10.1029/96JB00104 [URL last accessed 17/03/2025]



This document has been produced in the context of the Copernicus Climate Change Service (C3S).

The activities leading to these results have been contracted by the European Centre for Medium-Range Weather Forecasts, operator of C3S on behalf of the European Union (Delegation Agreement signed on 11/11/2014 and Contribution Agreement signed on 22/07/2021). All information in this document is provided "as is" and no guarantee or warranty is given that the information is fit for any particular purpose.

The users thereof use the information at their sole risk and liability. For the avoidance of all doubt , the European Commission and the European Centre for Medium - Range Weather Forecasts have no liability in respect of this document, which is merely representing the author's view.

Related articles


  • No labels