Conversion of time-varying Stokes coefficients into mass anomalies at the Earth’s surface considering the Earth’s oblateness

Time-varying Stokes coefficients estimated from GRACE satellite data are routinely converted into mass anomalies at the Earth’s surface with the expression proposed for that purpose by Wahr et al. (J Geophys Res 103(B12):30,205–30,229, 1998). However, the results obtained with it represent mass transport at the spherical surface of 6378 km radius. We show that the accuracy of such conversion may be insufficient, especially if the target area is located in a polar region and the signal-to-noise ratio is high. For instance, the peak values of mean linear trends in 2003–2015 estimated over Greenland and Amundsen Sea embayment of West Antarctica may be underestimated in this way by about 15%. As a solution, we propose an updated expression for the conversion of Stokes coefficients into mass anomalies. This expression is based on the assumptions that: (i) mass transport takes place at the reference ellipsoid and (ii) at each point of interest, the ellipsoidal surface is approximated by the sphere with a radius equal to the current radial distance from the Earth’s center (“locally spherical approximation”). The updated expression is nearly as simple as the traditionally used one but reduces the inaccuracies of the conversion procedure by an order of magnitude. In addition, we remind the reader that the conversion expressions are defined in spherical (geocentric) coordinates. We demonstrate that the difference between mass anomalies computed in spherical and ellipsoidal (geodetic) coordinates may not be negligible, so that a conversion of geodetic colatitudes into geocentric ones should not be omitted.


Introduction
Since the launch of Gravity Recovery and Climate Experiment (GRACE) satellite mission in 2002 (http://www.csr. utexas.edu/grace), satellite gravimetry has become one of the key tools to study large-scale mass transport in the Earth's system. Mass transport estimates based on GRACE data have been successfully applied in numerous studies of the solid Earth, cryosphere, oceans, and continental water resources (for an overview, see, e.g., Wouters et al. 2014). To produce these estimates, one exploits information about temporal gravity variations sensed by GRACE satellites.
Newton's law of gravitational attraction allows temporal gravity field variations to be computed uniquely, as soon as mass anomalies are given (see, e.g., Chao et al. 1987). Such a relationship is used, for instance, to esti-B Pavel Ditmar p.g.ditmar@tudelft.nl 1 Delft University of Technology, Stevinweg 1, 2628, CN, Delft, The Netherlands mate rapid gravity field variations caused by non-tidal mass transport in the ocean and atmosphere. These estimates are distributed as a so-called Atmosphere and Ocean De-aliasing product AOD1B , which is needed, along with other background force models, for a preprocessing of GRACE data. This is because rapid mass transport signals cannot be properly interpreted and play a role of additional source of noise in GRACE data. To properly process mass re-distribution in the atmosphere, it is essential to take into account its vertical structure (Swenson and Wahr 2002;Flechtner 2007), as well as the oblateness of the Earth (Forootan et al. 2013). All these considerations are taken into account in the production of the latest release (RL06) of the AOD1B product .
The inverse problem-conversion of temporal gravity field variations into mass anomalies-is more difficult. It cannot be solved uniquely without additional assumptions about the mass transport. A commonly used assumption is that mass transport takes place at a sphere of a certain radius a. Then, the Stokes coefficients describing the time-varying gravity field can be uniquely related to Fourier coefficients of mass anomalies, provided that the corresponding base functions are defined as surface spherical harmonics (Chao et al. 1987). Wahr et al. (1998) proposed to apply such an approach to recovery of mass anomalies from GRACE data, and this has been done since then in hundreds of studies.
Unfortunately, the overwhelming majority of authors ignore the fact that the actual Earth's surface substantially deviates from the sphere. If the radius a is defined as the mean equatorial radius or the semimajor axis of the reference ellipsoid (which is typically the case), these deviations become particularly large in polar areas. Obviously, such deviations may result in some distortions in mass anomaly estimates (Chao 2005(Chao , 2016. However, they are not quantified so far in the context of GRACE-based estimates, the temporal resolution of which is limited to about 400 km (in terms of wavelengths).
Furthermore, the commonly used expressions for a conversion of Stokes coefficients into mass anomalies refer to spherical (geocentric) coordinates of a point, including the geocentric colatitude θ . However, only a few authors define the exploited coordinates properly. In most of cases, θ is defined as just "the colatitude." We consider this as an indication that the step of converting standard ellipsoidal (geodetic) geographical coordinates into spherical ones is likely omitted in many studies, especially if they are conducted by nongeodesists. Though the difference between geocentric and geodetic colatitudes is relatively minor ( 11.5 ≈ 21 km), it might be improper to ignore it in all cases.
The primary goals of this publication are: (i) to demonstrate that mass anomaly estimates produced with the expressions from (Wahr et al. 1998) may not be sufficiently accurate (in particular, this concerns long-term mass losses in polar areas); (ii) to show that a simple modification of them may increase the accuracy by an order of magnitude; and (iii) to demonstrate that an (erroneous) interpretation of spherical coordinate θ as a geodetic colatitude may have a nonnegligible effect onto the estimated mass anomalies.
It is important to stress that our discussion is limited to a recovery of a 2-D mass re-distribution. In practice, this means that our focus is on hydrology, ice sheets, and oceans. The corresponding mass variations take place at the Earth's surface or just below it (at a depth not exceeding in most cases a few hundreds of meters). Then, our assumption that mass transport is confined to a thin layer near the Earth's surface is fully justified, and a quantification of it in terms of equivalent water heights is physically meaningful.
Mass transport in the solid Earth definitely occurs deeper, spanning a much larger range of depths. For instance, hypocenters of large earthquakes are typically located at the depth of a few tens of km, whereas glacial isostatic adjustment (GIA) takes place in the asthenosphere, the top boundary of which is located in most places at the depth of 100-200 km. Therefore, the techniques discussed in this paper are not applicable to mass transport in the solid Earth. A more extended discussion of limitations associated with a recovery of 3-D mass transport processes from GRACE data can be found in Chao (2016).
The structure of this paper is as follows. In Sect. 2, we present a general expression for the computation of mass anomalies at the Earth's surface, as well as possible simplifications of that expression, depending on the assumption about the Earth's surface shape. That section concludes by a comparison of the proposed simplifications using real GRACE data, which allows us to identify the expression that is the most appropriate in practice. In Sect. 3, we address the issue of the proper definition of the coordinate θ . Section 4 is left for a discussion and conclusions.

General information
Information delivered by GRACE is usually provided to the Earth science community in the form of monthly sets of Stokes coefficients. Those coefficients represent the mean value of the Earth's gravitation potential U (r , θ, φ, t) within a given month. The subtraction of a long-term mean value from each coefficient results in a time series of its temporal variations. Those variations can be linked to temporal variations U (r , θ, φ) of the gravitational potential (see, e.g., Heiskanen and Moritz 1967, Eqs. 2-39): where (r , θ, φ) are spherical coordinates of a given point (radial distance, colatitude, and longitude) in the terrestrial reference frame; G is the universal gravitational constant; M E is the Earth's mass; a is the Earth's radius (a more specific definition of this parameter is addressed below); l and m are the spherical harmonic degree and order, respectively; L max is a model-specific maximum degree; C lm and S lm are temporal variations of Stokes coefficients (by definition, S lm = 0 for m = 0); andP l,m are normalized associated Legendre functions. The argument of time is omitted in Eq. (1) for the sake of brevity. The summation does not contain the degree-0 term, since variations in the total mass of the Earth (which are described by that term) are negligible.
The exact definition of the Earth's radius a is a matter of convention. Classically, this parameter is defined as the semimajor axis of the reference ellipsoid or the mean Earth's equatorial radius: a ∼ 6,378,136 m (Heiskanen and Moritz 1967). GRACE data products are also provided in line with this definition. On the other hand, many publications define a, explicitly or implicitly, as the mean radius of the entire Earth: a = 6371 km (e.g., Chao et al. 1987;Wahr et al. 1998;Swenson and Wahr 2002). Strictly speaking, this means that an application of the expressions derived in those publications requires the corresponding rescaling of the GRACE-based Stokes coefficients. To the best of the author's knowledge, however, this is never done in practice. Apparently, the impact of this rescaling is assumed to be minor. In any case, in all the derivations presented below, a is defined as the semimajor axis of the reference ellipsoid.
To make the further derivations simpler, we rewrite Eq. (1) in a more compact form: (2) In this notation, spherical harmonic order m runs from −l to l andȲ lm are 4π -normalized surface spherical harmonics: The temporal variations C lm of Stokes coefficients in Eq.
(2) are related to the traditionally considered ones as In the further derivations, we make use of the publication by Swenson and Wahr (2002) as a starting point. We also use the same notation, when possible.
The general expression that connects temporal variations of density ρ(r , θ, φ) at/inside the Earth with temporal variations of Stokes coefficients is: where integration covers the entire sphere: whereas I l (θ, φ) describes vertically integrated density variations: Equation (5) is virtually equivalent to Eq. (2) in Swenson and Wahr (2002), with the exception that we watch the difference between the mean Earth's radius and the equatorial one.
As it is explained in Sect. 1, we assume that mass transport takes place in a thin layer at the Earth surface. As such, it can be represented by a single mass layer, so that I l (θ, φ) can be approximated as [cf. Eq. (7) in Swenson and Wahr (2002)]. In Eq. (8), σ (θ, φ) are variations of surface density (i.e., mass variations per unit area) and r s (θ, φ) is the function describing the shape of the Earth's surface that can be represented with a high accuracy as with ξ(θ, φ) being the height of the geoid above the sphere of radius a, whereas h(θ, φ) is the orthometric height of the Earth's surface topography [cf. Eq. (5) in Swenson and Wahr (2002)]. In the expressions below, the arguments in the functions ξ(θ, φ) and h(θ, φ) will not be explicitly written for the sake of brevity.
Variations of surface density can be related to the mass anomalies in terms of equivalent water heights (EWH) H w (θ, φ) in terms of equivalent water heights (EWH) with a simple formula where ρ w is water density. (10) Equations (5) and (8) describe the link between the temporal variations of surface density and the temporal variations of Stokes coefficients under the assumption that the Earth is a rigid body. In practice, solid Earth reacts elastically to changes in the load on the Earth' surface (see, e.g., Boy and Chao 2005). Hence, actual variations in the gravitational potential comprise both the direct effect of mass transport and the elastic deformation of the Earth (deformation potential). In order to take this into account, additional scaling factors (1 + k l ) are introduced, where k l are load Love numbers (Wahr et al. 1998), so that Eq. (8) turns into The function σ (θ, φ), as any other function of coordinates (θ, φ), can be represented in terms of the spherical harmonic expansion. After retaining the spherical harmonic degrees consistently with Eq. (1), we have: where C l m are Fourier coefficients, which can be computed on the basis of the Stokes coefficients C lm . To that end, we insert Eqs. (12) and (11) into Eq. (5), which yields: By interchanging the order of the summation and integration, we readily obtain: with Equation (14) represents a system of linear equations with constant coefficients B l,m,l ,m , which form a square matrix. By solving this system, one can transform the Stokes coefficients into the coefficients C l m . The system can be simplified further under some assumptions about the geometry of the Earth's surface, as discussed below.

Spherical Earth approximation (radius = 6378 km)
Let us assume that the Earth is the sphere of radius a, i.e., r s (θ, φ) = a, ξ(θ, φ) = 0, and h(θ, φ) = 0. Then, in view of the fact that the surface spherical harmonics form an orthogonal set on a sphere, the matrix formed by coefficients B l,m,l ,m turns into the unit one: where δ i,k is the Kronecker delta. Furthermore, the Earth's mass can be related to its mean density ρ E : where a E is the mean Earth's radius. Then, one can insert Eqs. (17) and (16) into (14). If the difference between a 3 E and a 3 is neglected (which is of the order of 0.3%), this readily results in: Thus, the computation of the coefficients C lm reduces to a scaling of the Stokes coefficients. A combination of Eqs. (18), (12), and (10) yields the well-known expression proposed for GRACE data processing by Wahr et al. (1998):

Ellipsoidal Earth approximation
Let the Earth's surface geometry be an ellipsoid of rotation.
is the height of the ellipsoid above the sphere of radius a. The radial distance r (θ ) of the points at the ellipsoidal surface is given by: where f is the ellipsoid flattening (WGS84 value: f = 1/298.2572) and e is eccentricity (e 2 = 2 f − f 2 ). Then, Eq. (15) turns into Furthermore, the trigonometric functions sin mφ and cos mφ form an orthogonal set in the interval [0; 2π ]: In view of Eq. (6), this allows the expression for the elements B l,m,l ,m to be simplified to: Thus, the system of linear equations given by Eq. (14) becomes block-diagonal and can be solved with ease.

Spherical Earth approximation (arbitrary radius)
Finally, one can also approximate the Earth's surface with a sphere of an arbitrary radius: , with r being an arbitrary value, provided that r > −a. In that case, the expression for the elements B l,m,l ,m simplifies to: so that the link between the coefficients C lm and the Stokes coefficients is: The parameter r can be chosen for a given target region such that the difference between the actual radial distance of the points at the Earth's surface and the Earth's equatorial radius a is taken into account. One step further is a "locally spherical" approximation. With this, we mean that the parameter r can be chosen not as a single constant for the entire target region, but as a value dependent on the location of the current point where the mass anomaly is computed. Let us assume, for example, that mass transport takes place at the ellipsoid, i.e., r = r (θ ) = ζ(θ). Then, a combination of Eqs. (25), (20), (12), and (10) yields the following expression for the mass anomalies in terms of EWH: or, in view of Eq. (17), This expression resembles Eq. (19), but contains an additional degree-depended scaling factor The most appropriate choice of the Earth's geometry approximation is discussed below.

Selection of the most appropriate approximation of the Earth's surface geometry
In the previous section, we presented a number of alternative expressions to convert temporal variations of Stokes coefficients into mass anomalies. The complexity of the associated computations depends on the adopted approximation of the Earth's shape. In order to identify the most appropriate computational scheme, we consider the computation of mass anomalies from real GRACE data.

Data
In this study, we use monthly gravity field solutions in 2003-2015 produced at the Center for Space Research (University of Texas at Austin) (Bettadpur 2012). For a few months, the solutions are absent. Furthermore, we ignored the solutions that were not limited to a specific calendar month. As such, 135 monthly solutions in total are exploited. Each of these solutions is formed by a set of Stokes coefficients complete to degree 96. No filtering is applied. It is expected that mass anomalies at relatively short spatial scales are particularly sensitive to the deviations of the assumed Earth's surface geometry from the actual one. Therefore, it makes sense to focus on the scenarios where the high-frequency signal in GRACE data is strong, whereas the noise level is low. In view of this, we use the time series of Stokes coefficients to estimate the mean rate of linear mass change (i.e., linear trend) in the entire time interval 2003-2015. In this way, random noise in the coefficients is largely suppressed. The linear trend is co-estimated with a bias and a quadratic term, as well as with annual and semiannual (co-)sinusoidal variations [see, e.g., Eq. (15) in Siemes et al. (2013), for more detail].

Computation of mass anomalies at the actual Earth's surface
As a reference, we convert the temporal variations of Stokes coefficients into mass anomalies using Eqs. (14) and (15) explicitly. To define the elevations h(θ, φ) of the Earth's surface, we use the Global Land One-kilometer Base Elevation (GLOBE) digital elevation model (GLOBE Task Team et al. 1999). For inland areas, the model provides terrain elevations above the mean sea level. The oceans are flagged; we set the elevations there equal to zero. The geoid heights ξ(θ, φ) above the sphere of radius a are approximated by Eq. (20) for the ellipsoid (thus, we ignore the differences between the reference ellipsoid and geoid as relatively minor, 100m). In general, the computed mass change rates are heavily contaminated by random noise, which is not surprising in view of the absence of filtering (not shown). Nevertheless, noise in polar areas is relatively low and allows one to clearly see strong signals over the territory of Greenland and the Amundsen Sea embayment of West Antarctica (Fig. 1). The strongest negative trends are observed at the Jakobshavn Isbrae (West Greenland) and at the Pine Island glacier (West Antarctica), see Fig. 1 and Table 1. A rapid ice mass loss at both locations is also revealed by other observation techniques (e.g., Groh et al. 2014;Mouginot et al. 2014). We explain the low noise level in Fig. 1, in spite of the absence of any filtering, by a combination of several factors. First, using 13 years of data allows random noise to be largely averaged out, as it is already mentioned above. Second, a large density of GRACE ground tracks, as well as the intersection of ascending and descending tracks at relatively large angles, ensures a good coverage of polar areas, which reduces random noise further. Third, Greenland and West Antarctica are notorious for the presence of very strong negative trends due to a rapid ice mass loss. Ironically, these locations are far away from the equator, so that deviations of the actual Earth's surface from the sphere of radius a ≈ 6378 km are large there. Thus, the presented areas can be considered as the "worst-case scenarios" for Eq. (19), which is traditionally used to convert Stokes coefficients into mass anomalies. Therefore, our further analysis is limited to the two geographical areas shown in Fig. 1.

Computation of mass anomalies at the sphere of 6378-km radius
Conversion of Stokes coefficients into mass anomalies with the commonly used Eq. (19)   from the values computed at the actual Earth's surface (see Fig. 2a). The observed differences are clearly anti-correlated with the total signal shown in Fig. 1. In other words, the traditionally used expression damps the recovered signals. This is not surprising, since the polar Earth's radius is 21 km smaller than the equatorial one, so that actual mass transport in the considered case takes place about 20 km further away from GRACE satellites than it is implicitly assumed in the traditional computations. The largest differences reach about 10 cm/year (Fig. 2a, Table 1). In order to better quantify the observed signal damping, we also represent the difference D(θ, φ) between the two variants of mass change rates in percentages, for which purpose the following expression is used:  is the mass change rate computed with the approximate formula andḢ (ref) w is the reference mass change rate. The results are presented in Fig. 2b. We only show the observed differences if the reference signal at a given point exceeds in absolute value 10 cm/year. In this way, we mask out the large relative differences that are caused by a small value in the denominator in Eq. (28). In spite of that, the observed relative differences are, in general, quite large. At the locations of strongest negative trends, they are of the order of 15% (Table 1), whereas at some other locations they are even larger, reaching 20% and more (Fig. 2b). One may argue that the largest relative differences may be associated with noise (e.g., over the ocean and, perhaps, over the inner part of Greenland). Mass anomaly estimates there must be filtered anyway, so that the errors introduced by the spherical Earth approximation are not critical. Nevertheless, even if we limit the discussion only to the coastal areas, the signal damping   Fig. 2, but the compared mass change rates are computed at the reference ellipsoid and at the actual Earth's surface caused by the considered approximation is at the level of up to 10-15%, which is hardly acceptable.

Computation of mass anomalies at the reference ellipsoid
Next, we compute mass anomalies under the assumption that mass transport takes place at the reference ellipsoid. To that end, we invert the block-diagonal matrix given by Eq. (23).
The resulting mass change rates are an order of magnitude closer to those computed at the actual Earth's surface, as compared to those produced under the assumption that the Earth is a sphere of 6378 km radius. The differences do not exceed 1 cm/year (Fig. 3a, Table 1). The relative differences are within 1.5% at the locations of the peak signal and typically stay within the 3% limit elsewhere (Table 1, Fig. 3b). Thus, the approximation of the Earth surface geometry with the reference ellipsoid may improve the conversion accuracy by an order of magnitude, as compared to the traditional approximation with the sphere of radius a = 6378 km. In addition, it is worth noticing that the spatial pattern of the observed differences shows a clear positive correlation with the signal itself (cf. Figs. 1 and 3a). In other words, the recovered signal is sharper than the actual one. Obviously, this is due to the fact that the surface of the reference ellipsoid (which is close to the sea level) is a few km further away from the GRACE satellites than the actual surface of ice sheets. Finally, relatively large differences are observed in the ocean areas at the western and southeastern coasts of Greenland. We relate them to a strong gradient of mass anomalies at the Greenland coasts. In view of a limited spectrum of the function H w (θ, φ), ringing artefacts associated with the Gibbs phenomenon must occur in those areas. As soon as recovered signal becomes sharper, these artefacts become more pronounced.

Computation of mass anomalies at the reference ellipsoid using the locally spherical approximation
Finally, we compute mass anomalies at the reference ellipsoid using the locally spherical approximation. That is, the spherical Earth's surface expression is used in the computations, but the radius of the sphere is latitude-dependent. At each latitude, it is set equal to the distance between the reference ellipsoid and the center of the Earth, cf. Eq. (26). The obtained linear trend estimates are surprisingly close to those produced with the explicit procedure addressed in the previous section. The differences between the results do not exceed 0.4 cm/year; see Fig. 4a. The relative differences stay at the level of at most 1-2% (Fig. 4b). Thus, the locally

Ellipsoidal versus spherical coordinates
In this section, we discuss the difference between mass anomalies computed in spherical (geocentric) and ellipsoidal (geodetic) coordinates. The maximum difference between geodetic and geocentric (co-)latitudes is observed near the 45 • latitudes, reaching approximately 11.5 arc-minutes or 21 km. In the polar areas, which are considered in our examples, this difference is smaller: of the order of 10-15 km or even less. One may argue that such differences must be negligible in the estimation of mass anomalies from GRACE data, since the spatial resolution of those estimates is a few hundreds of km. To demonstrate that such a statement may be unfair, we consider the meridional profiles that cross the Jakobshavn Isbrae in Greenland and the Pine Island glacier in West Antarctic (Fig. 5). The linear trend estimates along these profiles are presented as functions of both geodetic and geocentric latitude (Fig. 5). This figure clearly shows that the accuracy of locating a sharp signal may by far exceed the spatial resolution of GRACE data. As such, the difference between the estimates in spherical and ellipsoidal coordinates is clearly visible. Thus, we believe that the conversion of geodetic colatitudes into geocentric ones must not be omitted in GRACE data processing.

Discussion and conclusions
To convert time-varying Stokes coefficients into mass anomalies at the Earth's surface, geoscientists routinely use Eq. (19) or its equivalents. However, the results obtained with this expression represent mass transport at the spherical surface of 6378 km radius. In this study, we show that the accuracy of such a conversion may be insufficient, especially if the target area is located in a polar region and the signal-to-noise ratio is high. For instance, the mean linear trends in 2003-2015 estimated over Greenland and Amundsen Sea embayment of West Antarctica may be underestimated in this way by 10-15% or even more. Such an error may definitely exceed the noise level of current mass transport estimates. Moreover,  Fig. 1. The computed values are presented as functions of geocentric latitude (in red) and geodetic latitude (in green). In both cases, it is assumed that mass transport takes place at the actual Earth's surface, which is described by the GLOBE digital elevation model there are no doubts that the accuracy and spatial resolution of future estimates will increase further, so that a limited accuracy of the conversion based on Eq. (19) will likely become in the future even less tolerable. There are several reasons to expect improvements in the mass anomaly estimation in the foreseeable future: (i) an ongoing progress in the techniques for satellite gravity data processing; (ii) a continuously increasing duration of mass anomaly time series, which leads to a high accuracy of mean estimates for the entire available time interval (including mean linear trends); and (iii) the forthcoming launch of the GRACE Follow-On (GFO) satellite mission (https://gracefo.jpl.nasa.gov). GFO satellites will be equipped with a laser interferometer, allowing for an order-of-magnitude increase in the accuracy of inter-satellite ranging. Though the increase in the ranging accuracy may not result in a proportional increase in the accuracy of estimated mass anomalies in general , the level-1 data will definitely become cleaner at high frequencies, since ranging noise at those frequencies is dominant (Flury et al. 2008;Ditmar et al. 2012). Then, this reduction in noise level will likely have a positive effect onto the estimates of mass anomalies at small spatial scales, which are particularly vulnerable if the Earth's geometry is defined inaccurately.
As a solution, we propose an updated expression for the conversion of Stokes coefficients into mass anomalies. This expression is based on the assumptions that: (i) mass transport takes place at the Earth's surface that is approximated by the reference ellipsoid; (ii) at each point of interest, the Earth's surface is further approximated by the sphere with a radius equal to the current radial distance from the Earth's center ("locally spherical approximation"). The updated expression is nearly as simple as Eq. (19), but allows the inaccuracies associated with the conversion procedure to be reduced by an order of magnitude.
In addition, we demonstrate that it is advisable to convert geodetic (co-)latitudes into geocentric ones, when mass anomalies are computed. This is in spite of the fact that the shifts caused by this conversion are an order of magnitude smaller than the spatial resolution of GRACE-based estimates.
In "Appendix A," we summarize the recommended expressions for the conversion of Stokes coefficients into mass anomalies. Unlike in the main text, we use there the most traditional notation for the spherical harmonic expansion (when the orders start from 0 and not from −l), in order to facilitate the usage of the proposed expressions.
We would like to stress that the proposed conversion formula is particularly beneficial in the presence of strong signals in the range of high degrees, as it was the case in the considered examples. The absence of such signals makes the results much less sensitive to the assumption about the surface where mass transport takes place. For instance, the truncation of mass anomaly estimates at spherical harmonic degree 60 reduces the conversion errors 2-3 times, as compared to those presented in Fig. 2 (where the maximum degree is equal to 96).
In all the discussions so far, we assumed that the loading effects can be corrected just by introducing load Love numbers as additional scaling factors of the kind (1 + k l ). For a non-spherical Earth, such a simple approach is, strictly speaking, incorrect. However, the loading effects manifest themselves mostly in the range of low degrees. This can be understood from the fact that the scaling factors (1 + k l ) rapidly approach 1 as spherical harmonic degree increases. For instance, (1 + k l ) > 0.96 for degrees above 30 and (1+k l ) > 0.98 for degrees above 70 (Wahr et al. 1998). Since the impact of the proposed conversion formula is mostly limited to the range of relatively high degrees, we believe that a simplified treatment of the loading effects is justified.
Finally, one may pose the questions whether mass anomalies can be uniquely restored considering that mass transport takes place at the Earth's surface of a complicated geometry. For the case of a spherical Earth, a unique recovery of mass anomalies is guaranteed by Eq. (18). This formula establishes a unique link between the Stokes coefficients (that describe variations of gravitational potential) and the Fourier coef-ficients of mass anomalies. Moreover, this formula offers a practical way to make the corresponding conversion. By increasing the maximum spherical harmonic degree under consideration, one may, in principle, recover mass anomalies with an arbitrarily high spatial resolution. Thus, the spatial resolution of the results is fully defined by the spatial resolution of the exploited gravity field model. The situation with mass transport at the actual Earth's surface is more complicated. First of all, the unique conversion of time-varying Stokes coefficients into mass anomalies can only be guaranteed if the matrix composed of coefficients B l,m,l ,m (cf. Eq. 15) is invertible. In all the computations conducted in this study, this was indeed the case: this matrix was not only invertible, but also close to the unit one. However, it remains unclear if (or under what conditions) this matrix remains invertible in general. Furthermore, the presence of nonzero off-diagonal elements in this matrix implies that there is no unique link anymore between the spatial resolution of gravity field model and that of mass anomalies. For instance, high-frequency signals in terms of mass anomalies can map onto low-frequency signals in gravity field observations. If the former signals contain spherical degrees above L max in a given data processing run, a proper conversion of Stokes coefficients into mass anomalies becomes impossible. In other words, a realistic Earth's geometry may result in a new type of high-frequency signal aliasing, which is absent when mass transport takes place at a spherical surface. A quantification of this effect and, if necessary, designing optimal schemes to mitigate it are the subjects of further studies.