# A mass conserving filter based on diffusion for gravity recovery and climate experiment (GRACE) spherical harmonics solutions

Olivier Goux, J. Pfeffer, A. Blazquez, Anthony T. Weaver, M. Ablain

► To cite this version:

Olivier Goux, J. Pfeffer, A. Blazquez, Anthony T. Weaver, M. Ablain. A mass conserving filter based on diffusion for gravity recovery and climate experiment (GRACE) spherical harmonics solutions. Geophysical Journal International, 2023, 10.1093/gji/ggad016. hal-04720425

HAL Id: hal-04720425

https://hal.science/hal-04720425v1

Submitted on 4 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A mass conserving filter based on diffusion for gravity recovery and climate experiment (GRACE) spherical harmonics solutions

O. Goux<sup>1,3</sup>, J. Pfeffer,<sup>1</sup> A. Blazquez,<sup>2</sup> A. T. Weaver<sup>3</sup> and M. Ablain<sup>1</sup>

Accepted 2023 January 13. Received 2022 December 7; in original form 2022 March 15

## SUMMARY

Over the past two decades, the GRACE (Gravity Recovery and Climate Experiment) and GRACE Follow-On mission (GRACE-FO) have provided monthly measurements of the gravity field as sets of Stokes coefficients, referred to as spherical harmonics solutions. The variations of the gravity field can be used to infer mass variations on the surface of the Earth, mostly driven by the redistribution of water. However, unconstrained GRACE and GRACE-FO solutions are affected by strong correlated errors, easily identified as stripes along the north–south direction in the spatial domain. Here, we develop a filter based on the principle of diffusion to remove correlated errors and access the underlying geophysical signals. In contrast to many filters developed for this task, diffusion filters allow a spatially variable level of filtering that can be adapted to match spatially variable signal-to-noise ratios. Most importantly, the formalism of diffusion allows the implementation of boundary conditions, which can be used to prevent any flux through the coastlines during the filtering step. As mass conservation is enforced in the filter, global indicators such as trends in the global mean ocean mass are preserved. Compared with traditional filters, diffusion filters ensure the consistency of the solution at global and regional scales for ocean applications. Because leakage errors occurring during the filtering step are suppressed, better agreement is found when comparing diffusion-filtered spherical harmonic solutions with mascon solutions and independent estimates based on altimetry and *in situ* data.

**Key words:** Satellite gravity; Sea level change; Time variable gravity; Numerical modelling.

## 1 INTRODUCTION

Launched in 2002 and decommissioned in 2017, the Gravity Recovery And Climate Experiment (GRACE) is a mission composed of twin satellites following each other on a quasi-polar low Earth orbit (Tapley et al. 2004). A precise measure of the intersatellite distance, coupled with GPS and on-board accelerometer data, is used to estimate the monthly variations of the gravity field with an unprecedented accuracy (Tapley et al. 2019). Since June 2018, the GRACE Follow-On (GRACE-FO) mission continues the measurements of the gravity field with increased performance (Landerer et al. 2020), providing valuable insight on the mass transport within the Earth system over nearly two decades (Chen et al. 2021). Many geophysical phenomena contribute to the variations of the gravity field, including the redistribution of water in the atmosphere (Long et al. 2014), ocean (Landerer et al. 2015), hydrosphere (Rodell et al. 2018) and cryosphere (Chen et al. 2006), the deformation of the solid Earth due to surface excitation such as earthquakes (Panet et al. 2004) and deeper processes involving for example the viscous response of the Earth’s mantle to past deglaciation (Steffen et al. 2008).

The GRACE and GRACE-FO data have now become essential to monitor the global water cycle in a changing climate (e.g. Humphrey et al. 2016; Rodell et al. 2018; Pfeffer et al. 2022), and quantify critical indicators of climate change such as the Global Mean Ocean Mass (GMOM) increase (WCRP 2018; Barnoud et al. 2021) due to accelerated ice-mass losses in Greenland, Antarctica (Velicogna et al. 2020) and other inland glaciers (Cirací 2020).

To sustain these applications, several operational centres deliver monthly gravity field solutions consisting of a suite of Stokes coefficients of the gravitational potential up to a certain degree and order (Bettadpur et al. 2018), distributed for example by the International Centre for Global Earth Models (ICGEM) at http://icgem.gfz-potsdam.de/series (Ince et al. 2019). These solutions are usually referred to as Spherical Harmonics (SH) solutions. The Stokes coefficients representing the gravity field (i.e. Level 2 data) can then be projected onto the reference ellipsoid and converted into gridded surface mass anomalies (i.e. Level 3 data) expressed as equivalent water heights (EWH) anomalies with an expected resolution around 400 km (Ditmar 2018). Because of systematic instrumental errors and imperfections in the background

<footnote n="1"><sup>1</sup>Earth Observation, Magellium, 31520 Ramonville-Saint-Agne, France. E-mail: olivier.a.goux@gmail.com</footnote>

<footnote n="2"><sup>2</sup>EMC2, LEGOS, University of Toulouse, CNES, CNRS, UPS, IRD, 31400 Toulouse, France</footnote>

<footnote n="3"><sup>3</sup>ALGO-COOP, CERFACS / CECI CNRS UMR 5318, Toulouse, France</footnote>
models (Ditmar et al. 2012), GRACE-derived SH solutions are dominated by noise at high degrees and orders, resulting in characteristic elongated features in the north-south direction, commonly known as ‘stripes’ [shown in the panel (a) of Fig. 1]. As the stripes conceal the geophysical signals of interest, a post-processing step is required to remove correlated errors from SH solutions.

The first and most straightforward methods derived for that purpose are Gaussian filters (Jekeli 1981; Wahr & Molenaar 1998). As they need to be applied with a very large radius to efficiently remove the stripes, Gaussian filters tend to damp excessively the underlying geophysical signals. Moreover, one of the drawbacks of these filters is that they do not conserve the mass over naturally defined regions, such as the ocean (Chambers & Willis 2009), and more generally they attenuate the signal. In ocean applications, attenuation is largely due to the leakage of ground-based mass (i.e. ice, terrestrial water storage) variations into coastal waters during the filtering step, leading to an underestimation of the trend in the GMOM (Supporting Information, Fig. 4).

Since then, a wide range of more advanced filters has been developed, to improve the removal of stripe errors. However, such advanced filters do not attempt to reduce the leakage and fix the mass conservation issues observed with the Gaussian filter (Supporting Information, Fig. 4). The most common are presented hereafter. A more exhaustive review can be found in Crowley & Huang (2020). Kusche (2007) treated the SH solutions as realizations of a stochastic process which could be characterized through the modelled covariance matrix of the geophysical signal and a synthetic covariance matrix of the error. They defined in this framework an ‘optimal’ anisotropic filter, usually referred to as the DDK filter, provided as a large set of pre-computed coefficients directly applied on the Stokes coefficients. Eight sets of coefficients [originally four in Kusche et al. (2009)], corresponding to eight levels of filtering (also called order), are available, with DDK 8 corresponding to the lowest level of filtering and DDK 1 to the highest. The DDK3 filter is widely used for ocean applications and removes a large part of the noise from GRACE SH solutions (Fig. 1b). While the DDK filters assume that statistical properties of the geophysical signal and the error are time invariant, Horvath et al. (2018) developed a time-variable decorrelation (VADER) filter (often called VDK1,..., VDK8). The VADER filter discriminates more effectively between stripes and the geophysical signal than DDK filters (Fig. 1c) but requires more *a priori* information (i.e. full error covariance matrices for the expected signal and error for each month). This *a priori* information is not necessarily available for all SH solutions, and increases significantly the size of the data set. Moreover, users might want to preserve the independence of their filtered solutions from data sets that could be used for the validation of the final product. Alternatively, Wouters & Schrama (2007) attempted to reduce systematic errors in GRACE SH solutions using empirical orthogonal functions to isolate signals of geophysical interest from the noise. Other statistical decomposition techniques, such as independent components analysis (Frappart et al. 2011) or multichannel singular spectrum analysis (Prevost et al. 2019) were also applied in conjunction with more classical filters (i.e. Gaussian or DDK) to allow a better separation of geophysical signals from the remaining sources of errors.

Mass concentration (mascon) functions were developed to increase the temporal and spatial resolution of time-variable gravity solutions (e.g. Rowlands et al. 2005) and mitigate the noise present in the solutions without filtering (e.g. Luthcke et al. 2006). In the past years, a few mascon solutions have been released, for example

by JPL (Jet Propulsion Laboratory; Watkins et al. 2015), CSR (Center for Space Research; Save et al. 2016), GSFC (Goddard Space Flight Center; Loomis et al. 2019a) or ANU (Australian National University; Allgeyer et al. 2022). In mascon solutions, the Earth's surface is divided into grid elements, such as tiles, spherical caps (i.e. mascons) or irregular shapes following coastlines (Tregoning et al. 2022). In these elements, the mass is estimated to fit the range-rate or range-acceleration measurements (Level 1B) corrected for a mean background model including various geophysical phenomena (i.e. tides, ocean-atmosphere mass variations etc.). If unconstrained, the mascon inversion results in extremely noisy surface mass estimates, largely contaminated by stripe-like errors, similar to unfiltered SH solutions (e.g. Watkins et al. 2015). The regularization of the inversion is constructed to optimize spatially-variable signal-to-noise ratios, using either *a priori* constraints on expected geophysical signals (e.g. Watkins et al. 2015) or an iterative inversion procedure aiming at minimizing the residuals of the regression until the solution converges (e.g. Save et al. 2016; Loomis et al. 2019a). This results in an increased spatial resolution and reduced leakage errors (e.g. Watkins et al. 2015; Loomis et al. 2019a). However, mascon solutions come with processing choices defined by the processing centre (e.g. *a priori* constraints, background models, regularization matrices) that cannot be revisited for specific needs and applications. There is a much larger number of available Level 2 SH solutions, with a range of post-processing methods available, which allows a better assessment of the uncertainties linked with each choice (e.g. Blazquez et al. 2018).

As an alternative, we present in Section 3.1 a filter for SH solutions based on a diffusion operator that provides an anisotropic and spatially variable level of filtering, while mathematically enforcing the conservation of mass across any user-defined region (here, the ocean). The diffusion filter draws upon the framework of correlation operators based on diffusion, which are commonly used in data assimilation (Weaver & Courtier 2001) to model matrix-vector products with correlation matrices. The diffusion filter is implemented using sparse linear algebra techniques and is applied directly to gridded data. It is competitive in terms of computational cost with filters applied on Stokes coefficients. Its formalism allows the specification of boundary conditions to prevent the leakage of land signals into the ocean, as shown in Section 3.4. Section 3.5 shows that the diffusion filter enforces mass conservation within the boundaries and, in particular, the conservation of the GMOM. The diffusion filter can be tuned to apply an adequate level of filtering for each area, while preserving the global mass of continents, oceans, or any other user-defined area of interest. This allows the user to study local variations on land, on ocean, or at a global scale with a single filtered product in a consistent manner. This is illustrated through numerical experiments in Section 4, using data sets presented in Section 2. In Section 5, the different filtered and unfiltered solutions are compared in terms of variance in the spectral domain. Finally, Section 6 focuses on the ocean signal. The GMOM is estimated from filtered GRACE solutions and from estimates obtained with independent observations (altimetry and *in situ* data) in order to compare the solutions at a global scale. An evaluation of the remaining high-frequency noise in each solution based on the signal in the open sea is also presented. Finally, Section 7 summarizes the recommendations for the utilization of the diffusive filter.

A Python implementation of the diffusion filter used to produce the results of this article is available in open access at https://github.com/ogoux/diffusion_filter.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Linear trends from 2005 to 2015 of equivalent water height (EWH) in mm estimated from (a) SH unfiltered, (b) SH with a DDK filter of order 3, (c) SH filtered with a VDK filter of order 3 and (d) mascon. Note the different palette scales.</figcaption></figure>

## 2 DATA SETS

Four data sets are used throughout the present study. Our core data set is the Level 2 CSR SH solution, which is compared to the Level 3 CSR mascon solution and to an independent geodetic estimate of the ocean mass, based on the difference between the Sea Surface Height (SSH) measured by altimetry and steric sea level estimated using *in situ* measurements of the seawater temperature and salinity. The CSR data are available from 2002 to 2017 for Level 2 and Level 3, but with a lower quality after 2015 due to the lack of data from an accelerometer shut down on GRACE-B to prolongate the battery life of this satellite. The Level 2B GFZ data are used instead of the CSR data for the solution processed with a VADER filter, as this filter is not available for CSR data. On the other hand, the spatial coverage of *in situ* temperature and salinity data, based on Argo floats, is only satisfactory starting from 2005. SSH data from altimetry are available during the whole GRACE and GRACE-FO period.

We restricted our study period from January 2005 to December 2015 to allow meaningful comparisons between all data sources, in particular in terms of their contribution to the sea level budget.

Indeed, after 2016, the sea level budget closure might contain some issues independent of GRACE and GRACE-FO, such as the drifts of some Argo salinity captors (Barnoud *et al.* 2021) and Jason-3 radiometer (Barnoud *et al.* 2022).

### 2.1 GRACE SH solution

The GRACE SH solutions are monthly representations of the gravity field under the form of sets of Stokes coefficients. Here, we use the Release 06 of the Level 2 solution distributed up to the degree 60 by CSR (Bettadpur *et al.* 2018), except for the solution filtered with a VADER filter, which is a Level 2 solution distributed by the GFZ (Dahle *et al.* 2019) up to the degree 96 (as the VADER fitler is not available for the CSR solution). To allow meaningful comparisons with CSR mascon solutions, the same background models are used for the Level 2 and Level 3 solutions. The geocenter motion (degree 1 coefficients) is estimated from Sun *et al.* (2016). The oblateness coefficients (C20) are estimated after the SLR analyses by Loomis *et al.* (2019b). The C30 coefficients were not corrected, because both accelerometers were functional during the January 2005–December 2015 time period. The solutions are corrected for glacial isostatic
adjustment, using the ICE6G-D model by Peltier et al. (2016). The ocean circulation is restored in the SH solutions using the GAB model from the Release 6 of AOD1B. The final Stokes coefficients were converted to surface mass anomalies (expressed in equivalent water heights) and projected onto the WGS84 ellipsoid following Ditmar (2018), using the l3py package of Kvas (2018).

## 2.2 GRACE mascon solution

The Release 6 of the CSR mascon solution was used for comparison. For consistency with the Level 2 products, the GAD was removed from the mascon CSR solution and the GAB was restored. To facilitate comparisons between various data sets, we removed the mean over the reference period (2005–2015) for each of them.

## 2.3 Altimetry-based sea-level and Argo-based steric sea level

Monthly ocean mass changes are also estimated as the difference between the sea level changes from satellite altimetry and the steric sea level changes from an ensemble of in situ temperature and salinity measurements, based on Argo (WCRP 2018). Sea level anomalies (SLA) are estimated using the daily gridded data provided by the Copernicus Climate Change Service (C3S) at 0.25° resolution (Legeais et al. 2021). The C3S SLA data are averaged to produce monthly gridded data at a 1° resolution. A detailed description of the original C3S product is available in the Copernicus Climate Data Store.<fnref n="1" />

Steric sea level changes are estimated from monthly temperature and salinity grids provided at 1° resolution by the SIO (SCRIPPS Institution of Oceanography; Roemmich & Gilson 2009) as described in Marti et al. (2022). Following Purkey & Johnson (2010), a linear trend of +0.113 mm yr⁻<fnref n="1" /> representing the contribution of the deep ocean (between 2000 and 4000 m deep) is added to the ARGO-based steric SLA.

# 3 DIFFUSION BASED FILTERS

## 3.1 Principle of the diffusion filter

The physical process of diffusion naturally acts as a low-pass filter, as it tends to smooth and homogenize in space the quantities it affects. Numerical diffusion models are thus good candidates for a low-pass filter. One can show that, under the right conditions, the solution of the diffusion equation is equivalent to the application of a convolution filter. These properties allow the implementation of a flexible and computationally efficient filter based on classical numerical integration schemes used for solving the diffusion equation. The link between convolution filters and diffusion operators is exploited in Weaver & Courtier (2001) to design diffusion-based correlation operators for data assimilation algorithms. For that application, normalization operators are required to ensure that the implicitly modelled kernel of the correlation operator is a (correlation) function with unit amplitude. In our application, we ignore this criterion to enforce mass conservation instead.

The diffusion filter is described by a diffusion equation shown in eq. (1) for the special case where the diffusion coefficient $κ$ is constant.

$$ \frac{\partial \eta(t, z)}{\partial t} - \kappa \Delta \eta(t, z) = 0, \qquad (1) $$

$$ \eta(t = 0, z) = S(z), \qquad (2) $$

where $Δ$ denotes the Laplacian operator and $η$ denotes an arbitrary scalar function of a pseudo-time coordinate $t$ and a 2-D vector of spatial coordinates $z$. The quantity $S(z)$ denotes the signal to be filtered, and is provided as an ‘initial’ condition to the diffusion equation. The principle of diffusion filters is to integrate numerically eq. (1) up to a pseudo-time $T$ and define the filtered signal $S_f(z)$ as:

$$ S_f(z) = \eta(t = T, z). \qquad (3) $$

Both $S(z)$ and $S_f(z)$ do not depend on the coordinate $t$, which is referred to as a ‘time’ coordinate solely by analogy with diffusion. The filter is applied to a single set of gridded values at a given instant. In the current context, $t$ is related to the filter parameters as explained later. To compute $η(t = T, z)$, a numerical scheme has to be applied to eq. (1). In Weaver & Courtier (2001), an explicit Euler scheme is used. Explicit Euler schemes are straightforward to implement but are only stable if the discretization in time and space satisfies the Courant–Friedrichs–Lewy (CFL) condition. The CFL condition is restrictive as, for a fixed integration time $T$, it often imposes a large number of time steps. Instead, an implicit Euler scheme is generally preferable (Mirouze & Weaver 2010), because of its property of unconditional stability. In addition, this choice makes the filter more versatile, as the number of time steps acts as a degree of freedom to control the smoothing properties of the filter. In this respect, implicit diffusion is preferable for the filter.

Let $δ_t$ denote the length of a time step. The time-discrete counterpart of $η$ is

$$ \eta^k(z) = \eta(k\delta_t, z), \qquad (4) $$

where $k$ is a non-negative integer and $t = k\delta_t$. With an implicit Euler scheme, eq. (1) becomes

$$ \eta^{k+1}(z) = \mathcal{L}[\eta^k](z), \qquad (5) $$

$$ \mathcal{L}^{-1} \equiv \mathcal{I} - \delta_t \kappa \Delta, \qquad (6) $$

where $\mathcal{I}$ denotes the identity operator. The operator $\mathcal{L}$ represents one iteration of the implicit scheme (i.e. the propagation over one time step). Let $M$ be the number of time steps required to reach $T$, assuming that $δ_t$ is chosen to divide exactly $T$ (i.e. $T = Mδ_t$). Eqs (2), (3), (5) and (6), can then be combined:

$$ S_f(z) = \eta^M(z) = \mathcal{L}^M[\eta^0](z) = \mathcal{L}^M[S](z). \qquad (7) $$

The parameter $M$ and the product of $κ$ and $δ_t$ influence the nature and amount of filtering applied to $S(z)$. A property of interest of diffusion filters is that they are equivalent (when far enough from boundaries) to a convolution with functions from the Matérn class. These functions, commonly used to characterize stochastic processes (Whittle 1963; Guttorp 2006), vary in shape between a decaying exponential and a Gaussian function, as shown in Fig. 2.

The shape of the Matérn functions associated with diffusion filters is linked to the number of iterations $M$ and the product $κδ_t$, which can be interpreted as the square of a length scale. The Daley length scale $D$ is more intuitive to parametrize diffusion filters than the length scale $\sqrt{κδ_t}$. $D$ represents the half width at mid amplitude of the parabola osculating the Matérn function at its origin. $D$ defines the spatial extension of the function (Figs 2a and b). In two dimensions,

<footnote n="1">¹https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-sea-level-glo bal?tab=overview</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Influence of *D* and *M* on Matérn functions. (a) Effect of a variable Daley length scale *D* on a Matérn function for *M* = 4. (b) Fourier transforms of the Matérn functions shown in (a). (c) Effect of a variable number of iterations *M* on a Matérn function for *D* = 25 km. A Gaussian function with a radius of 25 km is shown for reference. (d) Fourier transforms of the Matérn functions shown in (c). All functions have been normalized to have an integral equal to one.</figcaption></figure>

*D* is linked to diffusion parameters through:

$$D = \sqrt{\kappa \delta_t (2M - 4)} \quad (8)$$

*M* is an integer that regulates the shape of the Matérn function (Figs 2c and d). When *M* tends to +∞ (with *D* fixed), it converges to a Gaussian function of radius *D*. As a consequence, if *D* is constant (and not spatially variable, anisotropic, or set to 0 locally as in the next sections), a diffusion filter would have an effect similar to a Gaussian filter of radius *D* for a large *M* (in practice *M* ≈ 8 – 10 is sufficient). *D* and *M* are sufficient to define the Matérn function. In practice, to use the definition of *D* from eq. (8), we need *M* to be greater than 2. For a given value of *D*, Matérn functions with different values of *M* have a comparable spatial extension (Fig. 2c). However, as *M* increases, they show a slower decay at the origin and a sharper decay at long range (Fig. 2c). While the effect of *M* is not obvious in the spatial domain, it appears clearly in the spectral domain (Fig. 2d). The Fourier transforms of the functions highlight their role as low-pass filters (Figs 2b and d). As a first approximation, we can consider that, while *D* determine the cut-off frequency of the filter (Fig 2b.), *M* determines the filter roll-off, that is how fast the signal is damped after this cut-off (Fig 2d.), and can thus be associated to the order of the low-pass filter.

In practice, the filter is applied to a field that is discrete in space. The spatial discretization and numerical implementation of the diffusion filter are presented in Appendix B.

## 3.2 Anisotropy

The suitability of the diffusion filter for the reduction of systematic errors in GRACE SH solutions is improved by using a more general form of the diffusion equation in which the constant diffusion coefficient is replaced by a spatially variable diffusion tensor *κ*; that is a symmetric, positive definite 2 × 2 matrix. The tensor can account for anisotropic filtering. It is introduced directly in the

Laplacian operator, which ensures that the diffusion smoothing kernel is symmetric (Weaver and Courtier 2001; Weaver and Mirouze 2013). Typically, a stronger level of filtering is needed in the east-west direction to remove the stripes in GRACE SH solutions while avoiding an unnecessarily strong damping of gradients of the geophysical signal along the north-south direction (Han et al. 2005). In a geographical coordinate system, this can be achieved using a diagonal diffusion tensor where different values are specified for the tensor elements associated with the zonal and meridional derivatives. The diagonal elements of *κ* correspond to the square of the Daley length scales in each direction, scaled by the factor 1/(2*M* – 4) that appears in eq. (8). This result in anisotropic Matérn functions with a larger spatial extension in the east-west direction than in the north-south direction, thus making it suitable for GRACE applications (Appendix B).

## 3.3 Spatially variable filtering

The diffusion filter may also be improved by introducing a spatially variable level of filtering. Typically, less filtering is needed over the continents than over the oceans, due to a stronger signal-to-noise ratio. To achieve a spatially variable level of filtering, the Daley length scale *D* can take different values can be applied to each grid cell. The filter is thus provided with a global field of *D*, from which spatially variable values of the diffusivity are deduced, with possibly a different value on each grid cell (to allow a spatially variable level of filtering) and each grid edge (to introduce anisotropy).

The possibility of applying spatially variable levels of filtering to SH solutions improves their consistency with mascon solutions. Indeed, the regularization used for mascon solutions allows a spatial variability in intensity. The DDK and VDK filters used for SH solutions are also spatially variable, but as they are defined in the spectral domain (i.e. in the spherical harmonics basis), it is more difficult to control their spatial filtering properties. In this study, the
parametrization was kept relatively simple, with only four different D values, varying for continents and oceans, as well as in north-south and west-east directions. A dependency on latitude has also been implemented for these four length scales on one of the filters, to account for the higher data coverage at high latitudes due to the polar orbit of the GRACE and GRACE-FO missions. The Daley length scales are modified by a factor that depends on the cosine of the latitude φ, and a user-defined scalar α as shown in eq. (9):

$$D_{\text{variable}}(\phi) = D_{\text{fixed}}(1 + \alpha \cos \phi). \qquad (9)$$

The Daley length scale $D_{fixed}$ specified by the user thus represents the minimal length scale reached at the poles. The factor α is used to set the magnitude of the latitude variability: if set to 0, the effective Daley length scale is equal to $D_{fixed}$ everywhere. At the Equator, the effective Daley length scale is inflated by a factor $1 + \alpha$.

## 3.4 Boundary conditions

To complete the filter definition, boundary conditions are needed. For geophysical applications on the surface of the Earth, periodic boundary conditions are used for the western and eastern boundaries of the grid. Grid cells have a trapezoidal shape, but grid cells adjacent to the Northern or Southern boundaries have one edge of zero length each, which makes the cells triangular. As the flux between two grid cells is proportional to the length of the edge separating them (like a 'surface' of exchange), no flux can go through the Northern and Southern boundaries. Boundary conditions preventing flux are usually referred to as Neumann boundary conditions, or reflective boundary conditions. Both periodic and Neumann boundary conditions are compatible with mass conservation on the domain: for western and eastern boundaries, any outbound flux on one side translates into an inbound flux on the other side; for the Northern and Southern boundaries, there is no flux at all through the boundaries.

While the previous boundaries are sufficient for the implementation of the operator, additional conditions can be specified to represent the physical properties of the field of interest. An example of interest for our application is the addition of coastlines. Indeed, a common issue in filtered GRACE solutions is the attenuation of the signal, in particular the leakage of hydrological signal from the continents into coastal waters. If we use a priori information on the positions of the coastlines (e.g. a land-sea mask), we can prevent any flux through the coastlines by modifying the diffusivity field, through the parametrization of the Daley length scale D. During the spatial discretization of the diffusivity tensor (Appendix A), a diffusivity value is attributed to each edge separating grid cells (Sections 3.2 and 3.3), deduced from D and M (eq. 8). Setting D to zero on an edge of a grid cell prevents any flux through this edge. Doing so on the whole coastline makes the land and ocean completely independent from each other. As a consequence, variations between the land and ocean signal are not damped by the filter. While coastlines are one of the most obvious example of such natural boundaries, any other delimitations may be specified depending on the application: for example, limits of drainage basins for hydrology. Similarly, Neumann boundary conditions can be used to avoid processing some areas. By setting D to zero on all edges around and inside the area, the input data in this area will stay unchanged and will not interact with its surroundings. We use Neumann boundary conditions to exclude areas without relevant measurements to compute the GMOM time series in Fig. 6. It can also be used to exclude

areas with invalid data so that it does not contaminate neighbouring valid data.

## 3.5 Conservation of the global mean ocean mass

This study aims at designing a filter that is able to reduce systematic errors in GRACE measurements by avoiding leakage of geophysical signals from the continental areas into coastal waters. In particular, we want to be able to interpret both regional ocean mass changes and global mean ocean mass changes. To preserve the GMOM through the filtering step, we need to define the limit between the ocean and the continents. The most straightforward way to do so is through a binary land/ocean mask. However, as our grid cells are relatively large ($1^{\circ} \times 1^{\circ}$), many of them contain significant surfaces of both land and water. On that account, a more accurate approximation can be obtained by using the ocean surface ratio of each grid each grid cell computed from the ETOPO1 global relief model. Let $\rho_i$ denote the ocean surface ratio of the ith element:

$$\text{GMOM}(\eta) = \frac{1}{S_o} \sum_i \eta_i s_i \rho_i, \qquad (10)$$

where $\eta_i$ denotes the global mass anomalies in ith grid cell, $s_i$ denotes the surface of the corresponding grid cell, and $S_o = \sum_i s_i \rho_i$ is the surface of the ocean. We can derive a condition for a general linear operator in matrix form A to preserve the GMOM:

$$\text{GMOM}(A\eta) = \text{GMOM}(\eta). \qquad (11)$$

Let $(A)_{i,j}$ denotes the element on the ith line and jth column of A. The left-hand side of eq. (11) can be rewritten as:

$$\sum_i \left( \sum_k (A)_{i,k} \eta_k \right) s_i \rho_i = \sum_k \left( \sum_i (A)_{i,k} s_i \rho_i \right) \eta_k. \qquad (12)$$

This term must be equated to the right hand side of eq. (11), which can be rewritten as:

$$\sum_k \eta_k s_k \rho_k = \sum_k (s_k \rho_k) \eta_k. \qquad (13)$$

As eq. (11) must hold for any vector $\eta$ (i.e. the GMOM must be conserved for an arbitrary input), we can obtain a simple criterion for A to conserve the GMOM:

$$\sum_i (A)_{i,k} s_i \rho_i = s_k \rho_k. \qquad (14)$$

The weighted average of the elements in the kth column of the operator A should then be equal to the weight corresponding to the kth grid cell.

Showing that the inverse of one step of diffusion conserves mass is sufficient to prove that the diffusion filter conserves mass. The inverse of one step of diffusion is a product with $WL^{-1}$ followed by a product with $W^{-1}$, and is thus equivalent to a product with $L^{-1}$ (see the discretization of the diffusion filter in Appendix A). One can show that the sum of the elements in the kth column of $WL^{-1}$ is equal to $w_k$ (this result can be deduced from eq. (76) of Mirouze & Weaver 2010, for example):

$$\forall k, \sum_i (WL^{-1})_{i,k} = w_k. \qquad (15)$$

The product with $W^{-1}$ is equivalent to dividing the elements of the ith line by $w_i$ so, from eq. (15):

$$\forall k, \sum_i (L^{-1})_{i,k} w_i = w_k. \qquad (16)$$
By comparing eqs (14) and (16), we see that the GMOM is conserved if the elements $w_i$ of W are equal to the product of the grid cell surface and ocean surface ratio $s_i \rho_i$. In order to accurately preserve the GMOM, mixed water/land cells should be weighted according to their ocean surface ratio.

This section provides a method to obtain a filtered solution that conserves the GMOM by mathematical construction. Consequently, diffusion filters conserve exactly the GMOM regardless of the values of $D$ or $M$ used, even in extreme cases: if $D$ is extremely large (e.g. larger than the perimeter of the Earth) on the ocean, the solution will progressively tend to a global average on the ocean (i.e. equal to the GMOM everywhere). On the other hand, if $D$ is close to zero, the solution will be not be modified and the GMOM will also be preserved. Here, the contribution of the ocean signal is assumed proportional to the ocean surface ratio. Continental signals due to hydrology or ice-mass changes are usually larger than the ocean signal, which may introduce errors in this approximation. To mitigate the influence of continental signals in mixed water/land cells, the ocean may be defined over a restricted area, including only with grid cells containing more than 90 per cent of ocean surface. This effectively reduces leakage, defined as the attribution of surface mass changes originated on land (hydrology, ice-melt) to the ocean. It is illustrated in Appendix C with examples of filtered solutions with two positions of the coastline (all mixed grid cells in the ocean; grid cells with more than 90 per cent of ocean surface in the ocean).

# 4 SURFACE MASS ANOMALIES FILTERED WITH DIFFUSION

In this section, we consider two different parametrizations of the diffusion filter referred to as the diffusion solutions A and B (Table 1). Each parametrization includes 4 length scale values for $D$, to allow different levels of filtering in the north-south and east-west directions over the continents and over the ocean. The solution A does not include any dependence on latitude to show how different regions are affected by the same length scale. The parameters of solution A have been chosen such that most of the noise is suppressed everywhere (at the expense of the geophysical signal in some regions). The parameters of solution B are latitude-dependent as described in eq. (9). The length scales provided for solution B in Table 1 are the length scales defined at the poles, and which are inflated by a factor $(1 + \alpha)$ at the Equator. The parameters of the diffusion filter of solution B have been chosen so that it provides the same level of filtering as the solution A around 30 degrees of latitudes, while using at higher-latitudes a level of filtering low enough so as not to damp excessively the geophysical signal (in particular in Greenland). It offers a good compromise between noise removal and preservation of the geophysical signal, and can be used as a starting point for most applications. The two diffusion solutions are compared to the same CSR SH solution filtered with a DDK3 filter, the CSR mascon solution and the GFZ solution filtered with a VDK3 filter (this filter is not available for CSR solutions).

The diffusion filter efficiently removes systematic errors identified as stripes in the unfiltered solutions (Figs 3a and b). Large scale geophysical signals, linked to the melting of ice-sheets (negative trends in Greenland and West Antarctica) and glaciers (negative trends in Alaska, Andes, Arctic Islands, Svalbard etc.), increasing the ocean mass (positive trends in most parts of the global ocean) are adequately recovered by both diffusion solutions, and are overall consistent with the mascon, DDK3 and VDK3 solution (Fig. 3). The trends recovered at large spatial scales in terrestrial water storage

are also consistent between the four solutions, with positive trends in the Amazon and Zambezi basins, as well as in East Australia, and negative trends over large parts of the United States of America, São Francisco river basin and Middle East regions (Fig. 3).

Compared to other SH solutions, the diffusion solutions A and B limit the leakage of mass signals from the continents into coastal waters, especially around strong ice melting points such as Greenland, the Antarctic peninsula, or the Gulf of Alaska. Large negative trends are observed in these coastal waters in the DDK3 and VDK3 solutions, which are reduced in the diffusion solutions A and B (Fig. 3). While diffusion filters prevent any leakage during the filtering step, some negative trends remain in the diffusion solutions (Figs 3a and b). Some of these negative trends, also visible in mascon solutions (and contributing to the negative trends of DDK/VDK solutions), are not necessarily linked to leakage and could be explained by the geophysical contribution of the sea level fingerprints (Tamisiea & Mitrovica 2011). The negative trends in the diffusion solutions could also be induced by a residual leakage contribution as diffusion filters completely prevent leakage during the filtering step but are not designed to reduce any previous contamination of the ocean signal. This contamination could be linked to the limited spatial resolution of the mass signals due to the ill-conditioned nature of the inverse problem linking the GRACE measurements to the mass signals. Negative trends of lower amplitude are also visible in the diffusion solutions in the Ross Sea near the Antarctic peninsula (Figs 3a and b). Part of these trends (just west of the Antarctic peninsula) are not visible in the mascon solution and could be caused by a contamination of the ocean signal by ice melt signal before filtering. Further west offshore, the negative trends observed are also visible in the mascon solution (Fig. 3e), which may indicate the emergence of significant geophysical signals in areas where large errors are expected due to leakage (Chen et al. 2022).

Globally, the spatial features observed on land in the diffusion solutions (Figs 3a and b) are consistent with those observed in the mascon solution (Fig. 3c). Because of the implementation of boundary conditions, both diffusion solutions preserve sharp transitions between the land and ocean, that are smoothed in the DDK3 or VDK3 solutions (see for instance the coasts of Northern Africa in Figs 4a and b, or the west coast of the United States in Figs 4f and g). The level of filtering applied homogeneously in solution A has been chosen to mitigate noise globally, and is thus appropriate for regions of low signal-to-noise ratio such as the Sahara, but excessive elsewhere, and especially in regions with a high signal-to-noise ratio such as Greenland. In this region, solution A fails to recover positive trends away from melting glaciers also visible in the DDK3, VDK3 and mascon solutions (as would DDK/VDK filters of order 1 or 2, as shown by Fig. S1). On the other hand, the positive trends inland Greenland are captured by solution B. The filter of solution B is latitude-dependent, and can thus been tuned to use the same level of filtering as solution A at the latitude of the Sahara, and a level of filtering low enough to preserve the geophysical signal at the latitude of Greenland. Smaller length scale values are indeed more adapted to high latitudes benefiting from a better coverage due to the polar orbit of the GRACE mission and large signal-to-noise ratios due to strong ice mass changes over continental areas. On the other hand, larger length scales values are clearly more adapted to tropical deserts such as the Sahara (Fig. 4). This discrepancy shows that there is no single level of filtering appropriate everywhere, and justifies the need for the latitude dependency used for solution B (and DDK/VDK solutions), which is able to damp most of the noise at low latitudes (Fig. 4) as well as preserve geophysical signal at high latitudes.
<table><caption>**Table 1.** Parameters of the diffusion solutions used in Figs 5 and 7. The subscript NS refers to the north-south direction while EW to the east-west direction. The length scales $D_{\text{ocean}}^{\text{NS}}$ and $D_{\text{ocean}}^{\text{EW}}$ are used on the ocean; the length scales $D_{\text{NS}}^{\text{land}}$ and $D_{\text{EW}}^{\text{land}}$ are used on land. The scalars $\alpha_{\text{NS}}$ and $\alpha_{\text{EW}}$ indicate, respectively, the extent to which extent length scales in the north-south and east-west directions are adjusted with the latitude as specified in eq. (9).</caption><thead><tr><th>Solution</th><th><i>D</i><sub>NS</sub><sup>ocean</sup> (km)</th><th><i>D</i><sub>EW</sub><sup>ocean</sup> (km)</th><th><i>D</i><sub>NS</sub><sup>land</sup> (km)</th><th><i>D</i><sub>EW</sub><sup>land</sup> (km)</th><th><i>M</i></th><th>α<sub>NS</sub></th><th>α<sub>EW</sub></th></tr></thead><tbody><tr><td>A</td><td>360</td><td>480</td><td>90</td><td>285</td><td>8</td><td>0</td><td>0</td></tr><tr><td>B</td><td>60</td><td>80</td><td>15</td><td>50</td><td>8</td><td>5.8</td><td>5.8</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Linear trends in surface mass anomalies from 2005 to 2015 estimated with the CSR SH solution filtered with the diffusion parametrization A (a) and B (b). The diffusion solutions A and B are compared to the CSR SH solution filtered with a classical DDK3 filter (c), the GFZ SH solution filtered with a VDK3 filter (d), and the CSR mascon solution (e). The parameters of the diffusion filters used for solutions A and B are given in Table 1.</figcaption></figure>

Overall, SH solutions filtered with diffusion seem more consistent with mascon solutions than the SH solutions filtered with DDK/VDK filters. Indeed, sharp variations of the trends through the coastlines are visible and the impact of hydrological and ice melt signals on the coastal waters is reduced. The difference of length scales between the land and the ocean allows a level of

filtering appropriate for both areas, without having to compromise between a noisy ocean signal and an excessively damped land signal. The diffusion filter used on solution B has been tuned to provide a good compromise between noise removal and preservation of the geophysical signal, and could be used as a recommended setting for a large range of applications requiring to reduce systematic errors
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Regional zoom of the linear trends of surface mass anomalies from monthly CSR SH and mascon solutions from 2005 to 2015. This figure uses the same parameters as Fig. 3.</figcaption></figure>

in GRACE solutions. However, the choice of optimal parameters is free and may change for specific needs. This versatility makes the diffusion filter adaptable to a large range of scientific applications, not limited to the reduction of systematic errors in GRACE data.

Beyond the choice of parameters such as *D* and *M* and $\alpha$, the position of the geographical boundaries of null mass fluxes has a strong influence on the filtered solution (see Appendix C). To mitigate the contamination of ocean signals by land signals, it is recommended to restrict the ocean area to grid cells where the ocean surface ratio is largely dominant.

## 5 SPECTRAL CONTENT OF THE SIGNAL

DDK, VDK and diffusion filters do not modify the power spectra of GRACE-based surface mass anomalies until degree 20. The power spectra is consistent with the mascon and unfiltered solutions (Fig. 5). Differences begin when the noise becomes predominant over the geophysical signal in the unconstrained solutions (Lenczuk et al. 2020). Around the degree 20, the amplitude of unfiltered solutions begins to grow exponentially as the noise progressively becomes predominant over the signal. The decay of power starting from degrees 30 to 40 in the diffusion solutions A and B confirms that the diffusion filter successfully reduces the high frequency noise. The power of the diffusion solution A and B decays slightly faster than the mascon solution (at degree 60, about 30 per cent lower for solution B and 40 per cent lower for solution A), which may be due to an attenuation of the geophysical signal. This is especially clear for the diffusion solution A which lacks small scale structures that are observed in the mascon solution and the diffusion solution B (see Section 4).

The power of the DDK 3 solution is close to the power of the mascon solution up to the degree 30 approximately, and drops quickly thereafter. It is much smaller than the other solutions at the highest degrees (about 95 per cent lower than the mascon solution at degree 60). The power spectrum of the VDK3 solution is larger but still lower than the other solutions at the highest degrees (about 80 per cent lower than the mascon solution at degree 60). On the other hand, the DDK3 and VDK3 solutions shows slightly more

power than the diffusion solutions at intermediate degrees (at most 15 per cent more than the solution B).

The removal of systematic errors identified as stripes proved to be efficient in the diffusion solution B and mascon solutions (Fig. 3), so that better preserved geophysical signals are likely to contribute to the relative gain of power at the highest degrees when compared to the DDK3 solution, including sharp transitions at the coastline (Fig. 5). However, the slightly larger power of DDK3/VDK3 solution at intermediate degrees could indicate that the diffusion solution B starts damping the signal too early (in terms of frequencies). One must keep in mind that the power spectrum is only an indicator of the average amplitude of spherical harmonic coefficients at a given degree. It does not allow us to differentiate the noise from the geophysical signal, and does not represent accurately the spatial and temporal variability of the gravity field, which suggests caution in the interpretation of the spectra.

## 6 OCEAN APPLICATIONS

The GMOM, which is computed as a global average over the ocean. In practice, to allow a comparison between all data sources, the ocean only includes areas defined such that:

* The ocean/land surface ratio is larger than 90 per cent.

* SLA data from satellite altimetry is available.

* Steric SLA from SIO is available.

To ensure that the comparisons are carried out on the same geographical domain, Neumann boundary conditions (described in Section 3.4) are used to exclude areas where altimetry or Argo data are not available (high latitudes, inner seas and a few coastal regions).

GMOM estimates computed from SH solutions filtered with diffusion filters match exactly the GMOM estimates computed from unfiltered SH solution by mathematical construction of the filter, regardless of which values are used for the Daley length scale *D* and the smoothness parameter *M*, as proved in Section 3.5. Because there is no attenuation over the oceanic domain due to the filter, the trend in the GMOM is exactly the same between the unfiltered solution and the diffusion solution (Fig. 6). GMOM changes are in-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Spectral analysis up to the degree 60 of the monthly GRACE solutions averaged over the period 2005–2015. The amplitude is shown in mm EWH. The parameters of the filters used for solutions A and B are summarized in Table 1.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Estimation of monthly Global Mean Ocean Mass (GMOM) changes from 2005 to 2015 using the CSR SH solution unfiltered and filtered with a Gaussian filter of radius 180 km (the geometric mean of the length scales of the diffusion filter A), the diffusion filter A and a DDK3 filter. Estimate of the GMOM are also evaluated from the GFZ SH solutions filtered with a VDK3 filter, the CSR mascon solution and the ocean mass derived from altimetry-based sea level and Argo-based steric sea level. Here the oceanic domain is defined with a restrictive geographical mask identical for all solutions considered.</figcaption></figure>

dependent of the set of parameters *D* and *M* of the filter. The GMOM retrieved with the diffusion solution is also found to be more consistent with the mascon solution and with the independent estimates computed with altimetry and ARGO data, than solutions filtered with a Gaussian, DDK or VDK filter (see Fig. 6 and Appendix D).

Indeed, a large reduction in the GMOM trend (–1 or –2 mm yr–1) is observed for solutions filtered with either a Gaussian filter or the more advanced DDK/VDK filters. Assessing whether the Gaussian filter reduces the trend more or less than the DDK/VDK solutions is difficult as there is no single filtering radius which would make it
<table><caption>**Table 2.** Recommended diffusion filter parameters for CSH SH solutions distributed up to the degree 60 and 96. The subscript NS refers to the north–south direction while EW to the east–west direction. The length scales $D_{\text{ocean}}^{\text{NS}}$ and $D_{\text{ocean}}^{\text{EW}}$ are used on the ocean; the length scales $D_{\text{NS}}^{\text{land}}$ and $D_{\text{EW}}^{\text{land}}$ are used on land. The scalars $\alpha_{\text{NS}}$ and $\alpha_{\text{EW}}$ indicate respectively to which extent length scales in the north–south and east–west directions are adjusted with the latitude as specified in eq. (9).</caption><thead><tr><th>Maximum degree</th><th>$D_{\text{NS}}^{\text{ocean}}$ (km)</th><th>$D_{\text{EW}}^{\text{ocean}}$ (km)</th><th>$D_{\text{NS}}^{\text{land}}$ (km)</th><th>$D_{\text{EW}}^{\text{land}}$ (km)</th><th>$M$</th><th>$\alpha_{\text{NS}}$</th><th>$\alpha_{\text{EW}}$</th></tr></thead><tbody><tr><td>60</td><td>60</td><td>80</td><td>15</td><td>50</td><td>8</td><td>5.8</td><td>5.8</td></tr><tr><td>96</td><td>80</td><td>100</td><td>40</td><td>60</td><td>8</td><td>4</td><td>6</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Global mean RMS of detrended and deseasoned ocean mass anomalies for unfiltered SH, diffusion A and B, DDK3 and mascon CSR solutions, and a VDK3 GFZ solution. The set of parameters of the diffusion filters used on solutions A and B are summarized in Table 1.</figcaption></figure>

comparable in terms of noise suppression to anisotropic, spatially variable filters such as DDK/VDK or diffusion filters.

A common practice when estimating the GMOM is to mask coastal waters within a certain distance from land to limit the influence of the leakage (Chambers & Willis 2009). The application of a mask might in turn bias the GMOM and increase its uncertainty, as coastal waters are not equivalent to open seas in terms of redistribution of water mass (Tamisiea & Mitrovica 2011). Because the diffusion filter does not introduce leakage in the filtered solution, the gridded solutions can be used directly to compute the GMOM without needing a coastal mask.

Another useful indicator on the ocean is the RMS on the ‘quiet’ ocean, that is where polar areas have been excluded between -60° and +60° of latitude, as well as ocean areas within 1000 km from the coast, to remove known sources of geophysical signals (e.g. Circumpolar Antarctic Current, Arctic Oscillation, mass transport in shallow seas). The resulting mask is shown in Appendix E. Linear trends, annual cycles and semi-annual cycles are removed in this area using the outputs of a classical ordinary least squares inversion at each grid cell. Assuming that ocean mass changes are negligible at interannual timescales away from extra-tropical ocean basins and shallow seas, the RMS of the detrended and deseasoned ocean mass anomalies provides a good proxy of the level of noise in GRACE solutions. The RMS of unfiltered solutions ranges between 100 and 300 mm, indicating the presence of large systematic errors in GRACE solution (Fig. 7). With a DDK3 or VDK3 filter, the expected level of noise over ‘quiet’ areas of the ocean goes down to 20–30 mm. The expected level of noise is further reduced when using a mascon or a diffusion solution, and is generally comprised

within 10–20 mm, which indicates a very high performance in the removal of systematic errors in GRACE SH solutions. Increasing further the length scales in the parametrization of the diffusion filter would reduce the RMS of detrended and deseasoned ocean mass anomalies even further as the solution would progressively tend to a global average. However, the diffusion solutions and the mascon solution already reach a low level of RMS, where the contribution of the residual geophysical signal is likely significant. Indeed, while expected to be small in tropical ocean basins, water mass changes may occur at interannual time scales, reaching up to 20–25 mm for known climate modes such as the El Nino Southern Oscillation or North Atlantic Oscillation (Pfeffer et al. 2022). Below this level, reducing the RMS is thus not necessarily a good indicator of noise suppression.

With the diffusion filter, it is therefore possible to interpret both regional and global ocean mass changes in a consistent way. Here, we limited the spatial extent of the ocean to be consistent between all observation systems (satellite radar altimetry, satellite gravimetry and *in situ* T/S measurements), but the same conclusions can be reached with other coastline positions (see Appendix D for the GMOM trends with different coastline positions) or filter parameters ($D$ and $M$ values).

## 7 RECOMMENDATIONS

In this paper, we used the CSR solution distributed up to the degree 60 to illustrate the effect of diffusion filters, and provided recommended filter parameters (Solution B in Table 1). Diffusion filters can also be applied to higher degree solutions, such as the CSR solution distributed up to the degree 96, though the filtering radii would need to be adapted. As high degree coefficients contain more noise, larger filtering radii are required to efficiently damp the high frequency noise (Table 2). Regional trends obtained with the degree-96 solution are shown in the supporting information, showing that diffusion filters can successfully remove noise in the degree-96 solution, while preserving the GMOM. However, increasing the filtering radius also limits the added value of higher degree coefficients. The parameters of the diffusion filter recommended in Table 2 offer a fair compromise between noise removal and preservation of the underlying geophysical signal, and can be used for most applications. Beyond the classical corrections summarized in Section 2.1 (e.g. GAB, C20, etc.), no pre-processing is required before applying a diffusion filter. Pre-filtering the solution with a Gaussian or DDK filter [as done for other post-processing methods, e.g. Frappart et al. (2011) or Prevost et al. (2019)] should be avoided to prevent the introduction of leakage errors. A Gaussian filter can be seen as a particular case of diffusion filters. With a constant Daley length scale $D$ in both directions, a large value of $M$ ($\ge 8$) and no boundary conditions enforced within the domain (i.e. on coastlines), the effect of a diffusion filter becomes equivalent to that of a Gaussian filter of radius $D$. Diffusion filters can thus act as a pre-filtering method before using another post-processing method instead of a Gaussian
filter, in order to enforce mass conservation, introduce anisotropy and/or spatial variability of the filtering radius.

## 8 CONCLUSION

Diffusion filters provide a new post-processing method to reduce noise in SH solutions while allowing for the conservation of mass across any user-defined geographical domain, as well as spatially variable levels of filtering that can be adapted to different signal-to-noise ratios. We demonstrated here that the properties of diffusion filters can be adapted to avoid leakage through the coastlines and conserve exactly the global mass of the ocean independently of the level of filtering applied. This allows for a globally consistent estimate of surface mass anomalies with different levels of filtering across continents and oceans. Diffusion filters are computationally efficient and can be easily applied on any SH solution projected onto the spatial domain, without the need for prior information about the error statistics of SH solutions.

The diffusion filter is based on the numerical solution of a diffusion equation, with an effect similar to a convolution filter. In this study, the effective smoothing kernel of the diffusion filter provides stronger filtering in the east-west direction to target more effectively the stripe-like errors in the GRACE-derived SH solutions. The length scales in each direction are defined as spatially variable fields. The number of diffusion iterations acts as an additional shape parameter that determines the order of the filter. At large orders, the effect of the diffusion filter would be close to that of a Gaussian filter if the user-specified length scales were isotropic and spatially constant, but other shapes of convolution kernels can be obtained with lower orders. A Python implementation of the diffusion filter has been made available in open access at https://github.com/ogoux/diffusion_filter.

While the definition of boundary conditions prevents leakage through the coastline during the filtering step, it has no effect on leakage errors already present in the Level 2 solution, linked to the limited spatial resolution of the GRACE measurements. Furthermore, many grid cells contain both land and ocean surface, which can introduce some land signal in the ocean (or vice versa) depending on the position of the coastline. For a better separation of land and ocean signal, alternative land masks could be considered.

As this article mainly serves as a proof of concept, we have kept the filter parameter specifications rather simple: two length scales are defined for the land, and two for the ocean (one per direction). These length scales are inflated at low latitudes in one of the diffusion solutions to account for the difference of data coverage associated with the GRACE orbit at low and high latitudes. This dependency on latitude and the separation of land and ocean allows the level of filtering to match the variations of the signal-to-noise to first order.

This study was conducted using a SH solution distributed up to the degree 60, but diffusion filters can also be applied to solutions available up to higher degrees. For these solutions, a new set of parameters is required in order to remove the additional high-frequency noise at high degrees. We provided a set of parameters for SH solutions available up to the degree 96, which can be used to damp the high-frequency noise in these solutions while conserving the GMOM. More advanced parametrizations could help to extract high-resolution geophysical signals from higher degree coefficients affected by a stronger and more heterogeneous level of noise than lower degree coefficients. The diffusion filter parameter specifications could be improved by defining a set of length scales that match

more closely the expected geophysical signal-to-noise ratio beyond the dependency to latitude, for example by discriminating between biomes (e.g. lower filtering in ice covered areas, higher filtering in tropical deserts, etc.). These settings could draw from the methods used to tune the intensity of the regularization of mascon solutions. The parameters could also be used to approach a statistically optimal filter using the same approach as DDK or VDK filters. These advanced parametrizations would allow a better separation of geophysical signals from the noise, in particular at high degrees where this noise is prevalent.

The filter can also be made more versatile by using a linear combination of SH solutions filtered with different parameters: as long as the coefficients of the linear combination sum to one, the mass conservation of the signal is still enforced. Such combinations of filters are equivalent to a convolution filter whose kernel can show features such as negative lobes that are often observed in the kernels of statistical filters such as DDK (Appendix C). Alternatively, using different filtering length scales for each of the *M* iterations of the filter also opens up a wider range of kernels, whose spectrum show a steeper slope after the cut-off frequency.

Finally, diffusion filters are not limited to the applications addressed in this paper. Boundaries can be defined to delimit any area such as hydrological basins, glaciers or coastlines as in this study. Moreover, while the filter has been designed to reduce systematic errors in surface mass anomalies derived from GRACE SH solutions, the current implementation can be used to filter other gridded geophysical data sets, such as geoid or altimetry data sets. On those data sets, the filter could be parametrized to conserve other conservative quantities than the GMOM (by providing other weights *w* than the ocean/land surface ratio). Diffusion filters could also be extended to unstructured data sets using methods already in use for correlation operators (Guillet et al. 2019).

## ACKNOWLEDGMENTS

We would like to thank the editor and the reviewer for their remarks which proved helpful for improving this article. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (GRACEFUL Synergy Grant agreement No 855677).

## DATA AVAILABILITY

The GRACE monthly SH and mascon solutions are published by the CSR on the PODAAC drive: https://podaac-tools.csr.nasa.gov/drive/files/GeodeticsGravity/grace. The C3S daily altimetry grids are available on the Climate Data Store: https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-sea-level-global?tab=overview. The SIO ARGO temperature and salinity gridded products are available at: http://sio-argo.ucsd.edu/RG.Climatology.html.

## REFERENCES

* Allgeyer, S. et al., 2022. ANU GRACE data analysis: orbit modeling, regularization and inter-satellite range acceleration observations, *J. geophys. Res.*, **127**(2), e2021JB022489, doi:10.1029/2021JB022489.

* Barnoud, A. et al., 2021. Contributions of altimetry and Argo to non-closure of the global mean sea level budget since 2016, *Geophys. Res. Lett.*, **48**(14), e2021GL092824, doi:10.1029/2021GL092824.

* Barnoud, A., Pfeffer, J., Cazenave, A. & Ablain, M., 2022. Revisiting the global mean ocean mass budget over 2005–2020, *EGUsphere* [preprint], https://doi.org/10.5194/egusphere-2022-716, 2022.
Bettadpur, S., 2018. Level-2 Gravity Field Product User Handbook, CSR.

Blazquez, A., Meyssignac, B., Lemoine, J.M., Berthier, E., Ribes, A. & Cazenave, A., 2018. Exploring the uncertainty in GRACE estimates of the mass redistributions at the Earth surface: implications for the global water and sea level budgets, *Geophys. J. Int.*, **215**, 415–430.

Chambers, D.P. & Willis, J.K., 2009. Low-frequency exchange of mass between ocean basins, *J. geophys. Res.*, **114**(C11), doi:10.1029/2009JC005518.

Chen, J.L., Wilson, C.R. & Tapley, B.D., 2021. Satellite gravity measurements confirm accelerated melting of Greenland ice sheet, *Science*, **313**, 1958–1960.

Chen, J.L., Tapley, B., Tamisiea, M.E., Save, H., Wilson, C., Bettadpur, S. & Seo, K.W., 2021. Error assessment of GRACE and GRACE follow-on mass change, *J. geophys. Res.*, **126**(9), doi:10.1029/2021JB022124.

Chen, J.L., Cazenave, A., Dahle, C., Llovel, W., Panet, I., Pfeffer, J. & Moreira, L., 2022. Applications and challenges of GRACE and GRACE follow-on satellite gravimetry, *Surv. Geophys.*, **43**, 304–345.

Ciraci, E., Velicogna, I. & Swenson, S., 2020. Continuity of the mass loss of the world's glaciers and ice caps from the GRACE and GRACE follow-on missions, *Geophys. Res. Lett.*, **47**(9), doi:10.1029/2019GL086926.

Crowley, W.J. & Huang, J., 2020. A least-squares method for estimating the correlated error of GRACE models, *J. geophys. Res.*, **221**, 1736–1749.

Dahle, C. et al., 2019. The GFZ GRACE RL06 monthly gravity field time series: processing details and quality assessment, *Remote Sens.*, **11**, 2072–4292.

Ditmar, P., 2018. Conversion of time-varying Stokes coefficients into mass anomalies at the Earth's surface considering the Earth's oblateness, *J. Geod.*, **92**, 1401–1412.

Ditmar, P., Teixeira da Encarnação, J. & Hashemi Farahani, H., 2012. Understanding data noise in gravity field recovery on the basis of inter-satellite ranging measurements acquired by the satellite gravimetry mission GRACE, *J. Geod.*, **86**, 441–465.

Frappart, F., Ramillien, G., Leblanc, M., Tweed, S.O., Bonnet, M.P. & Maisongrande, P., 2011. An independent component analysis filtering approach for estimating continental hydrology in the GRACE gravity data, *Remote Sens. Environ.*, **115**(1), 187–204.

Guttorp, P. & Gneiting, T., 2006. Studies in the history of probability and statistics XLIX; on the Matérn correlation family, *Biometrika*, **93**, 989–995.

Guillet, O., Weaver, A.T., Vasseur, X., Michel, Y., Gratton, S. & Gürol, S., 2019. Modelling spatially correlated observation errors in variational data assimilation using a diffusion operator on an unstructured mesh, *Q. J. R. Meteorol. Soc.*, **145**, 1947–1967.

Horvath, A., Murböck, M., Pail, R. & Horwath, M., 2018. Decorrelation of GRACE time variable gravity field solutions using full covariance information, *Geosciences*, **8**(9), doi:10.3390/geosciences8090323.

Humphrey, V., Gudmundsson, L. & Seneviratne, S.I., 2016. Assessing global water storage variability from GRACE: trends, seasonal cycle, subseasonal anomalies and extremes, *Surv. Geophys.*, **37**, 357–395.

Ince, E.S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F. & Schuh, H., 2019. ICGEM – 15 years of successful collection and distribution of global gravitational models, associated services, and future plans, *J. Geophys. Res.*, **114**(2), 647–674.

Jekeli, C., 1981. Alternative Methods to Smooth the Earth's Gravity Field, Dep. of Geod. Sci. and Surv., Technical Report No. 327, Geodetic and GeoInformation Science, Department of Civil and Environmental Engineering and Geodetic Science, The Ohio State University Columbus, OH.

Kusche, J., 2007. Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models, *J. Geod.*, **81**, 733–749.

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.

Kvas, A., 2018. I3py package, available at https://doi.org/10.5281/zenodo.1450900.

Landerer, F.W., Wiese, D.N., Bentel, K., Boening, C. & Watkins, M.M., 2015. North Atlantic meridional overturning circulation variations from GRACE ocean bottom pressure anomalies, *Geophys. Res. Lett.*, **142**(19), 8114–8121.

Landerer, F.W. et al., 2020. Extending the global mass change data record: GRACE follow-on instrument and science data performance, *Geophys. Res. Lett.*, **12**, 647–662.

Legeais, J-F., Meyssignac, B., Faugère, Y., Guerou, A., Ablain, M., Pujol, M-I, Dufau, C & Dibarboure, G., 2021. Copernicus sea level space observations: a basis for assessing mitigation and developing adaptation strategies to sea level rise, *Front. Mar. Sci.*, **8**, 2296–7745.

Long, D., Longuevergne, L. & Scanlon, B.R., 2014. Uncertainty in evaporation from land surface modeling, remote sensing, and GRACE satellites, *Water Resour. Res.*, **50**(2), 1131–1151.

Loomis, B.D., Luthcke, S.B. & Sabaka, T.J., 2019a. Regularization and error characterization of GRACE mascons, *J. Geod.*, **83**, 1381–1398.

Loomis, B.D., Rachlin, K.E. & Luthcke, S.B., 2019b. Improved Earth oblateness rate reveals increased ice sheet losses and mass-driven sea level rise, *Geophys. Res. Lett.*, **46**, 6910–6917.

Lenczuk, A., Leszczuk, G., Klos, A. & Bogusz, J., 2020. Comparing variance of signal contained in the most recent GRACE solutions, *Geod. Cartogr.*, **69**, 69–37.

Luthcke, S.B., Rowlands, D.D., Lemoine, F.G., Klosko, S.M., Chinn, D. & McCarthy, J. J., 2020. Monthly spherical harmonic gravity field solutions determined from GRACE inter-satellite range-rate data alone, *Geophys. Res. Lett.*, **33**(2), doi:10.1029/2005GL024846.

Marti, F. et al., 2021. Monitoring the ocean heat content change and the Earth energy imbalance from space altimetry and space gravimetry, *Earth Syst. Sci. Data*, **14**, 229–249.

Mirouze, I. & Weaver, A.T., 2010. Representation of correlation functions in variational assimilation using an implicit diffusion operator, *Q. J. R. Meteorol. Soc.*, **136**, 1421–1443.

Panet, I. et al., 2004. Coseismic and post-seismic signatures of the Sumatra 2004 December and 2005 March earthquakes in GRACE satellite gravity, *Geophys. J. Int.*, **171**(1), 177–190.

Peltier, W.R., Argus, D.F. & Drummond, R., 2016. Comment on “An Assessment of the ICE-6G-C,(VM5a) GlacialIsostatic Adjustment Model” by Purcell et al., *J. geophys. Res.*, **123**, 2019–2028.

Pfeffer, J., Cazenave, A. & Barnoud, A., 2022. Analysis of the interannual variability in satellite gravity solutions: detection of climate modes fingerprints in water mass displacements across continents and oceans, *Clim Dyn*, **58**, 1065–1084.

Prevost, P., Chanard, K., Fleitout, L., Calais, E., Walwer, D., van Dam, T. & Ghil, M., 2019. Data-adaptive spatio-temporal filtering of GRACE data, *Geophys. Res. Lett.*, **219**(3), 2034–2055.

Purkey, S.G. & Johnson, G.C., 2010. Warming of Global Abyssal and Deep Southern Ocean Waters between the 1990s and 2000s: contributions to global heat and sea level rise budgets, *J. Clim.*, **23**, 6336–6351.

Rodell, M., Famiglietti, J.S., Wiese, D.N., Reager, J.T., Beaudoing, H.K., Landerer, F.W. & Lo, M.-H., 2018. Emerging trends in global freshwater availability, *Nature*, **557**, 651–659.

Roemmich, D. & Gilson, J., 2009. The 2004–2008 mean and annual cycle of temperature, salinity, and steric height in the global ocean from the Argo Program, *Prog. Oceanogr.*, **82**, 81–100.

Rowlands, D.D., Luthcke, S.B., Klosko, S.M., Lemoine, F.G.R., Chinn, D.S., McCarthy, J.J., Cox, C.M. & Anderson, O.B., 2005. Resolving mass flux at high spatial and temporal resolution using GRACE intersatellite measurements, *Geophys. Res. Lett.*, **32**(4), doi:10.1029/2004GL021908.

Save, H., Bettadpur, S. & Tapley, B.D., 2016. High-resolution CSR GRACE RL05 mascons, *J. geophys. Res.*, **121**(10), 7547–7569.

Han, S.C., Shum, C.K., Jekeli, C., Kuo, C.Y., Wilson, C. & Seo, K.W., 2005. Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement, *Geophys. J. Int.*, **163**(1), 18–25.

Sun, Y., Riva, R. & Ditmar, P., 2016. Optimizing estimates of annual variations and trends in geocenter motion and J2 from a combination of GRACE data and geophysical models, *J. geophys. Res.*, **121**(11), 8352–8370.

Steffen, H., Denker, H. & Müller, J., 2008. Glacial isostatic adjustment in Fennoscandia from GRACE data and comparison with geodynamical models, *J. Geodyn.*, **46**(3–5), 155–164.
Tamisiea, M.E. & Mitrovica, J.X., 2011. The moving boundaries of sea level change: understanding the origins of geographic variability, *Oceanography*, **24**, 24–39.

Tapley, B.D., Bettadpur, S., Ries, J.C., Thompson, P.F. & Watkins, M.M., 2004. GRACE measurements of mass variability in the Earth system, *Science*, **305**(5683), 503–5.

Tapley, B.D. et al., 2019. Contributions of GRACE to understanding climate change, *Nat. Clim. Chang.*, **9**(5), 358–369.

Tregoning, P., McGirr, R., Pfeffer, J., Purcell, A., McQueen, H., Allgeyer, S. & McClusky, S.C., 2022. ANU GRACE data analysis: characteristics and benefits of using irregularly shaped mascons, *J. geophys. Res.*, **127**(2), doi:10.1029/2021JB022412.

Velicogna, I. et al., 2020. Continuity of ice sheet mass loss in Greenland and Antarctica from the GRACE and GRACE follow-on missions, *Geophys. Res. Lett.*, **47**(8), doi:10.1029/2020GL087291.

Wahr, J. & Molenaar, M., 1998. Time variability of the Earth's gravity field: hydrological and oceanic effects and their possible detection using GRACE, *J. geophys. Res.*, **103**, 205–229.

Watkins, M.M., Wiese, D.N., Yuan, D., Boening, C. & Landerer, F.W., 2015. Improved methods for observing Earth's time variable mass distribution with GRACE using spherical cap mascons, *J. geophys. Res.*, **120**, 2648–2671.

WCRP Global Sea Level Budget Group, 2018. Global sea-level budget 1993-present, *Earth Syst. Sci. Data*, **10**, 1551–1590.

Weaver, A.T. & Courtier, P., 2001. Correlation modelling on the sphere using a generalized diffusion equation, *Q. J. R. Meteorol. Soc.*, **127**, 1815–1846.

Weaver, A.T. & Mirouze, I., 2013. On the diffusion equation and its application to isotropic and anisotropic correlation modelling in variational assimilation, *Q. J. R. Meteorol. Soc.*, **139**, 242–260.

Whittle, P., 1963. Stochastic processes in several dimensions, *Bull. Inst. Int. Stat.*, **40**, 974–994.

Wouters, B. & Schrama, E. J. O., 2007. Improved accuracy of GRACE gravity solutions through empirical orthogonal function filtering of spherical harmonics, *Geophys. Res. Lett.*, **34**(23), doi:10.1029/2007GL032098.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Linear trends from 2005 to 2015 of equivalent water height (EWH) in mm estimated from (a) CSR SH with a DDK filter of order 1, (b) CSR SH with a DDK filter of order 2, (c) GFZ SH with a VADER filter of order 1 and (c) GFZ SH with a VADER filter of order 2.

**Figure S2.** Linear trends from 2005 to 2015 of equivalent water height (EWH) in mm estimated from (a) the CSR SH solution up to the degree 96 with a DDK filter of order 1, (b) the CSR SH solution up to the degree 96 with a DDK filter of order 3 and (c) the GFZ SH solution up to the degree 96 with a VADER filter of order 3 and (c) the CSR mascon solution.

**Figure S3.** Linear trends from 2005 to 2015 of equivalent water height (EWH) in mm estimated from (a) the CSR SH solution up to the degree 96 with a DDK filter of order 1, (b) the CSR SH solution up to the degree 96 with a DDK filter of order 3 and (c) the GFZ SH solution up to the degree 96 with a VADER filter of order 3 and (c) the CSR mascon solution.

**Figure S4.** Linear trend of the GMOM over the period 2005–2015 estimated from different filtered SH solutions. The CSR SH solution distributed up to the degree 60 is filtered with a Gaussian filter, a diffusion filter and a DDK filter. The GFZ solution distributed up to the degree 96 is filtered with a VDK filter.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the

authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## APPENDIX A: SPATIAL DISCRETIZATION AND IMPLEMENTATION OF THE DIFFUSION FILTER

The discrete counterparts of the functions $S, S_f$ and $\eta^k$ are the vectors $S, S_f$ and $\eta^k$. Note that 2-D fields are represented by 1-D vectors. The matrix $L^{-1}$, defined as the spatially discretized counterpart of $L^{-1}$, is used to implement the operator described by eq. (7):

$$L^{-1} = I - \delta_t \nabla \cdot (\kappa \nabla), \quad (A1)$$

where $I$ denotes the identity matrix, $\nabla$ the divergence operator, $\nabla$ the gradient operator, $\kappa$ the spatially variable diffusivity tensor and $\Delta$ the Laplacian discretized here with finite differences in curvilinear coordinates (Weaver & Courtier 2001). Curvilinear coordinates offer a natural framework to account for the variations of the size and shape of the grid cells with latitude when the grid lies on a sphere or an ellipsoid.

The application of one 'time step' of the filter is then equivalent to solving a linear system:

$$L^{-1} \eta^{k+1} = \eta^k, \quad (A2)$$

where the matrix $L^{-1}$ is positive definite and self-adjoint with respect to a diagonal matrix $W$ whose elements are the surfaces of each grid cell. The self-adjoint matrix on the left hand side of (A.2) can be converted into a symmetric matrix by left-multiplying (A.2) by $W$. This makes (A.2) easier to solve using standard matrix solvers. The application of the filter can be summarized as in Algorithm 1:

**Algorithm 1** Diffusion filter

Compute $\kappa$ from the inputs $D$ and $M$ ($\delta_t$ is fixed to 1);

Build the non-zero elements of $WL^{-1}$ and $W$;

Initialize $\eta^0 = S$;

for $k$ from 1 to $M$ do

Compute the product $W\eta^{k-1}$;

Solve $(WL^{-1})\eta^k = W\eta^{k-1}$ for $\eta^k$;

end for

Fetch the filtered value $S_f = \eta^M$.

The low computational cost of this method stems from the sparsity of the matrices involved. Indeed, apart from a few elements displaced by boundary conditions, the matrix $WL^{-1}$ is penta-diagonal. For a grid of resolution 1°, less than 0.01 per cent of the elements of $WL^{-1}$ are different from zero. Consequently, if a sparse storage format is used to generate $WL^{-1}$, and a sparse solver for symmetric positive definite matrices is used for the system, the procedure becomes less demanding in terms of memory than building explicitly the matrix of a convolution product. The computational cost is approximately proportional to $M$, but remains low enough to be competitive with filters applied on spherical harmonics. As an order of magnitude, for $M=4$, about 190 s in CPU time are needed to filter a thousand 1° grids simultaneously (multiple grids can be filtered simultaneously to vectorize computations). As the code relies on the built-in parallelization routines of Scipy, the actual computation time is likely to be lower on most architectures (e.g. about 30 s of elapsed time on a four core laptop).
## APPENDIX B: INTRODUCTION OF ANISOTROPY IN THE DIFFUSION KERNEL

The convolution kernels associated with the diffusion operator are not specified explicitly. However, they can be visualized by applying the filter to a 'Dirac' field, that is a vector equal to zero everywhere except for one pixel where it is equal to 1. We can also set multiple pixels to 1 (far enough apart so as not to interact together) in order to compare two different convolution kernels: one in the ocean near the western shore of Ireland, one at the same latitude on land in central Europe (Fig. A1), and one 20° of latitude lower in the middle of Algeria. We can also use this method for filters applied on SH coefficients directly such as DDK or Gaussian filters by converting this Dirac field into SH coefficients, applying the filter, and reprojecting the result on the ellipsoid.

The diffusion filters used on the solution A and B (see Table 1 for the specific choice of parameters) both use anisotropic convolution kernels, with a stronger level of filtering in the east-west direction (Figs A1a and b). Both solutions have a weaker level of filtering on land than on the ocean, resulting in convolution kernels with a smaller spatial extent. The convolution kernels are restricted either to the ocean or land due to the boundary conditions enforced in the diffusion filter on coastlines.

On the other hand, filters applied to SH coefficients such as DDK3 or Gaussian filters have similar convolution kernels over land and ocean (see Figs A1c and d). A convolution kernel from these filters can cover both land and ocean at the same time, which results in a leakage through the coastlines. The convolution kernels of the diffusion filter B and the DDK 3 filter decrease in size with latitude. In practice, the distortion due to the projection, visible for the diffusion filter A or the Gaussian filter whose kernels appear different even though they are not latitude-dependent, can hide this property. It is also important to note that the kernels are projected on a grid of fixed resolution in terms of longitude and latitude. The amplitude kernel covering a given physical distance will thus be

'diluted' over more grid points in the east-west direction at higher latitudes.

## APPENDIX C: DEFINITION OF THE OCEAN/CONTINENT BOUNDARY

Beyond the choice of the parameters *D* and *M*, the position of the coastline also has an influence in the diffusion filter, because of the no flux condition imposed at the land/ocean boundary. Here, we compare two CSR SH solutions filtered with diffusion, using the same parameters *D* and *M* (see *D* and *M* values for solution B in Table 1), and differing only by the position of the coastline. In Fig. A2(a), only the grid cells where the ocean surface ratio is larger than 90 per cent are included in the ocean (as in the solution B used in previous sections). In Fig. A2(b), all grid cells containing some ocean surface are included in the ocean. In diffusion solutions, leakage occurs due geophysical signals originating from land displaced on the 'ocean side' of the coastline before the filtering step. As the solution of panel (b) includes more land in mixed land/ocean cells, an additional leakage is expected in this solution. It is visible for example off the shore of Greenland just north of Iceland, or around the Antarctic peninsula. The weighting by the ocean surface ratio of the mixed land/ocean grid cells mitigates this phenomenon, yet the difference in amplitude between the ice melt signal and the surrounding ocean response makes the remaining continental signal preponderant over the ocean signal. To avoid this phenomenon, a more restrictive definition of the ocean is needed as done in panel (a)

## APPENDIX D: TREND OF THE GMOM TIME SERIES

## APPENDIX E: AREA OF LOW VARIABILITY OF THE OCEAN
<figure><img src="image_10.png" /><figcaption>**Figure A1.** Examples of smoothing kernels of the diffusion filters used for solutions A and B (see Table 1) are shown in panels (a) and (b), respectively. The equivalent convolution kernels for a DDK3 filter and a SH Gaussian filter with a radius of 300 km are shown for comparison in panels (c) and (d), respectively.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure A2.** Linear trends of a CSR SH solution filtered with two different diffusion filters. Both filters use the parameters from the solution B in Table 1 but different positions for the coastlines. Grid cells with at least 90 per cent of ocean surface are included in the ocean in panel (a). All grid cells with some ocean surface are included in the ocean in panel (b).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure A3.** Mask used in Section 6 to define the ‘quiet’ ocean over which the spatial RMS is estimated and shown in Fig. 7.</figcaption></figure>

<table><caption>**Table D1.** Linear trends of the estimation of the monthly Global Mean Ocean Mass (GMOM) from 2005 to 2015 using the CSR SH solution unfiltered and filtered with Gaussian filters of radii 180 and 480 km, the diffusion filter A and a DDK3 filter. The GMOM is also estimated from the GFZ SH solutions filtered with a VDK3 filter, the CSR mascon solution and the ocean mass derived from altimetry-based sea level and Argo-based steric sea level. The two different radii of the Gaussian filter correspond, respectively, to the geometric mean of the length scales used in the solution A (180 km), and the largest length scale used in the solution A (480 km). Two different positions are considered for the coastline: one including all mixed land/ocean grid cells in the ocean, one including only grid cells with more than 90 per cent of ocean surface.</caption><thead><tr><th>Coastline position</th><th>0 per cent of<br>ocean surface</th><th>90% per cent<br>of ocean<br>surface<br>(mm yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>SH, unfiltered</td><td>1.811</td><td>1.811</td></tr><tr><td>SH, gaussian (180 km)</td><td>1.711</td><td>1.712</td></tr><tr><td>SH, gaussian (480 km)</td><td>1.502</td><td>1.500</td></tr><tr><td>SH, DDK3</td><td>1.671</td><td>1.674</td></tr><tr><td>SH, VDK3</td><td>1.611</td><td>1.615</td></tr><tr><td>mascon</td><td>1.847</td><td>1.846</td></tr><tr><td>Altimetry - ARGO</td><td>1.811</td><td>1.803</td></tr><tr><td>SH, diffusion</td><td>1.811</td><td>1.811</td></tr></tbody></table>