A review of different mascon approaches for regional gravity field modelling since 1968

The geodetic and geophysical literature shows an abundance of mascon approaches for modelling the gravity field of the Moon or Earth on global or regional scales. This article illustrates the differences and similarities between the methods, which are labelled as mascon approaches by their authors.

Point mass mascons and planar disc mascons were developed for modelling the lunar gravity field from Doppler tracking data. These early models had to consider restrictions in observation geometry, computational resources or geographical pre-knowledge, which influenced the implementation. Mascon approaches were later adapted and applied for the analysis of GRACE observations of the Earth's gravity field, with the most recent methods based on the simple layer potential.

Differences among the methods relate to the geometry of the mascon patches and to the implementation of the gradient and potential for field analysis and synthesis. Most mascon approaches provide a direct link between observation and mascon parameters – usually the surface density or the mass of an element – while some methods serve as a post-processing tool of spherical harmonic solutions. This article provides a historical overview of the different mascon approaches and sketches their properties from a theoretical perspective.

1 Introduction

The gravity field of the Earth influences daily life in many ways. Local plumb lines define the upward direction, and several scientific instruments must be levelled before usage. Gravity measurements provide corrections for geophysical height systems and allow the exploration of mineral deposits or caves. A gravity field model is also required for inertial navigation systems within aeroplanes, ships or submarines. On a regional or global scale, mass redistributions – like the melting of glaciers or the changes in ground water – are reflected in the temporal variations of the gravity field.

The gravity field of the Earth – or another celestial body – can be analysed on a global scale when enough orbiting satellites are tracked, even without special onboard instruments or ground-based measurements. Satellite data provide a fast and homogeneous sampling in contrast to ground-based observations. The gravity field analysis establishes a connection between the tracking data and a set of base functions to model the gravity field. An analysis by spherical harmonic functions is preferable for spherical bodies, as these base functions are the natural solution of the Laplace equation. This set of base functions is also a complete and orthogonal system, which simplifies the analysis. However, a reasonable spherical harmonic analysis requires orbit observations with global and homogeneous data distribution, and the model will have the same resolution everywhere.

Alternative localising base functions – e.g. point masses, spherical radial base functions, wavelets, Slepian function – are investigated and applied when the data distribution is irregular or when more details in a region of interest shall be detected. This article will summarise the localising base functions, which are labelled as mascons and applied for the gravity field modelling of the Earth and Moon.

Studies of the early lunar orbiters demonstrated significant orbit disturbances, which were traced back to an irregular lunar gravity field. The term “mascon” was introduced by Muller and Sjogren (1968b) for describing these mass concentrations near the surface. In the same work, the name mascon was also introduced for the mathematical modelling of these mass concentrations. The concept was applied for several years to the gravity field of the Moon, as the method was capable of the nearside restriction of data in contrast to spherical harmonic solutions. Interest in regional modelling of the Earth's gravity field has increased significantly since the gravity field mapping mission GRACE (2002–2017) and its successor mission GRACE-FO (2018–present). The new observations enabled the analysis of temporal variations caused by the redistribution of water and ice masses, where regional gravity field modelling overcomes the spherical harmonic solutions. Hence, the mascon concept has been adapted and applied to Earth-related data by several research groups, either for regions of interest (Luthcke et al., 2008; Schrama et al., 2014; Ran et al., 2018) or on a global scale (Koch and Witte, 1971; Andrews et al., 2015; Save et al., 2016).

A closer inspection of the publications, however, shows a variety of approaches under the label of mascons. This article will give a historical overview of the most prominent representatives and an adequate definition of the mascon base functions. All different meanings of the investigated mascon approaches can be covered by the following definition: the term mascon either refers to the fact of a significant gravitational anomaly within a celestial body or to a modelling of these anomalies by localising base functions. The localising base functions, which are labelled as mascons, include point masses or discrete surface elements based on the simple layer potential. In the case of surface elements, the surface density is constant per mascon, and each localising base function is – in a spectral representation, at least in the limit of high-degree expansion – a two-dimensional step function on the sphere. Methods of post processing are also labelled as a mascon approach when their surface elements have a constant surface density. The shape of the mascon is not relevant for the definition, and the surface of the celestial body is not necessarily covered.

This publication is focused on the mascons' definitions and will ignore other processing steps, like background models, regional constraints or regularisation techniques. Each mascon approach is presented by the associated gravitational potential of a single element and its gradient in the notation of representative literature. The properties of each approach are deduced from the theoretical perspective only, but without treating programming experiments or numerical aspects. Such a detailed and comprehensive review of the different mascon approaches cannot be found in literature to our knowledge.

In several previous articles, the authors quote only the original publication (Muller and Sjogren, 1968b) for the term mascon and restrict themselves in the following texts to a specific mascon approach with its literature (e.g. Luthcke et al., 2008; Lemoine et al., 2007; Krogh, 2011; Andrews et al., 2015).

A point mass model and planar discs are applied for modelling the lunar gravity field in Wong et al. (1971), and both methods are considered as mascon approaches. In Watkins et al. (2015) and Save et al. (2016), different mascon approaches are presented in the introductions, but without formulas or historical background. The authors of both articles classify three principal concepts:

• A.

  mascons that have an analytical expression for the gravitational potential and explicit partial derivatives for the gradient;

• B.

  mascons that are represented by a finite series of spherical harmonic functions and with partial derivatives derived via the chain rule;

• C.

  mascons that serve as a post-processing tool to obtain regional mass changes from monthly spherical harmonic solutions.

An analogous classification with additional literature is presented by Abedini et al. (2021a), whose contribution is a numerical method for the gradient, which does not fit into the threefold scheme.

Many recent publications are related to the mascon solutions of either the NASA Goddard Space Flight Center (GSFC), the Jet Propulsion Laboratory (JPL) or the Center for Space Research (CSR). The current JPL solutions are spherical cap mascons with analytical partial derivatives – i.e. category A in the classification – which are presented in Sect. 3.2. The mascon approaches of GSFC and CSR are based on spherical harmonic functions, and they are a prominent example of type B (see Sect. 3.1). The mascon visualisation tool at the University of Colorado Boulder (https://ccar.colorado.edu/grace/index.html, last access: 23 September 2022) enables an analysis and comparison of the latest solutions at JPL and GSFC for regions and generates time series per location.

2 Mascons for modelling the lunar gravity field

2.1 Mascons – mass anomalies close to the Moon's surface

The origin of the mascon concept is closely related to early models of the lunar gravity field.

In the space race between the Soviet Union (USSR) and the United States of America (USA), both nations wanted to send their representatives to the Moon first. The possible landing sites were investigated by spacecrafts, starting with Luna 1 (USSR) in 1959, which missed the Moon due to navigation issues. The first man-made object on the Moon was the space probe Luna 2 (USSR) in a design impact in 1959, followed by several missions by both nations. The spacecraft Luna 10 (USSR) and Lunar Orbiter 1 (USA) were the first artificial orbiters around the Moon in 1966 (Neal, 2008).

In both orbiter missions, the observed orbits differed after a short time from the predicted ones, which indicated either an incorrect or an incomplete model. As other error sources could be excluded soon, the orbit disturbances were explained by significant mass anomalies below the Moon's surface. For these anomalies, the term “mass concentration” or “mascon” was introduced in Muller and Sjogren (1968b).

Figure 1 Lunar maria and selected tracks of the Lunar Orbiter given in Muller (1972).

All identified mascons on the nearside of the Moon cause relatively large and positive effects up to 200 mGal, and their locations are one-to-one correlated to the major lunar maria, including Imbrium, Serenitatis, Crisium, Humorum and Nectaris, which are visualised in Fig. 1 (Muller, 1972).

In particular for the Moon it is still common to call a large area with a significant positive mass anomaly a mascon independent of its mathematical representation (Floberghagen, 2001, p. 3). A similar behaviour can be found, for example, in Barthelmes (1986, p. 35), where the mass anomalies of the Earth's gravity field are called mascons without using the phrase for the mathematical modelling as well. This thesis focuses on the point mass modelling, but it also sketches simple layer potential with discrete surface elements, and both aspects will be identified as mascon approaches in the current article.

2.2 Point mass mascons

A quick modelling of the mass anomalies was important for the preparation of the latter space missions and the landing on the Moon. The chosen representation should

• consider the geographical pre-knowledge, i.e. the lunar maria as expected locations of the mass anomalies;

• consider the observation geometry, i.e. the fact that only the near side of the Moon allows observations from terrestrial ground stations;

• enable a direct relation between observables – Doppler tracking data in the case of the early lunar missions – and the estimated mascon parameters;

• remain simple due to limited computer resources.

The first three requirements are still important arguments for regional gravity field analysis – by mascons, wavelets, radial basis functions, Slepian functions, etc. – while the limited resources implied a simple modelling of the anomalies by point masses.

The original papers (Muller and Sjogren, 1968a, b; Muller, 1972) lack a formula representation of the potential, but it is re-constructed, for example, in Floberghagen (2001, p. 19):

$$\begin{array}{}\text{(1)}& V\left({\mathit{r}}_{\mathrm{P}}\right)=GM\left({\displaystyle \frac{\mathrm{1}}{\Vert {\mathit{r}}_{\mathrm{P}}\Vert }}-\sum _{q=\mathrm{1}}^Q{\displaystyle \frac{\mathit{\delta }{m}_q}{\Vert {\mathit{r}}_{\mathrm{P}}-{\mathit{r}}_q\Vert }}\right),\end{array}$$

with

• V(r_P): gravitational potential at the calculation point r_P;

• G: gravitational constant;

• M: mass of the celestial body;

• δm_q: mass ratio between point masses and total mass M;

• r_q: centres of the point masses.

Please note that, for consistency, all mascon quantities and their geometries are labelled in this article by means of an index (here: $q=\mathrm{1},\mathrm{2},\mathrm{\dots },Q$), and the calculation point is labelled by the index P, both independent of the cited articles.

2.2.1 Relation to the observation and the estimation process

A standard observation technique for space probes is the Doppler tracking, i.e. the change in frequency of a (re)-transmitted signal due to the relative motion of the spacecraft and the ground station. The American missions use a few globally distributed stations, which meanwhile form the Deep Space Network of the NASA and which are operated by JPL today¹. The Doppler signal does not provide complete information on the position or velocity; rather, it only projects the relative velocity between station and space probe onto the line of sight (Muller and Sjogren, 1968a; Weinwurm, 2004; Floberghagen, 2001).

The relationship between observation and mascon parameters requires a description of the change in the velocity – i.e. the acceleration – of the spacecraft caused by the gravitational potential. Hence, it is sufficient to derive the gradient by ${\ddot{\mathit{r}}}_{\mathrm{P}}=\mathrm{\nabla }V\left({\mathit{r}}_{\mathrm{P}}\right)$ of the potential. For point mass models, the gradient is calculated via

$$\begin{array}{}\text{(2)}& \mathrm{\nabla }V\left({\mathit{r}}_{\mathrm{P}}\right)=GM\left(-{\displaystyle \frac{\mathrm{1}}{\Vert {\mathit{r}}_{\mathrm{P}}{\Vert }^{\mathrm{3}}}}{r}_{\mathrm{P}}+\sum _{q=\mathrm{1}}^Q{\displaystyle \frac{\mathit{\delta }{m}_q}{\Vert {\mathit{r}}_{\mathrm{P}}-{\mathit{r}}_q{\Vert }^{\mathrm{3}}}}\left({\mathit{r}}_{\mathrm{P}}-{\mathit{r}}_q\right)\right).\end{array}$$

To emphasise the special requirements and restrictions for the early lunar modelling, some details will be sketched here as well: according to Muller and Sjogren (1968a, b), residual observations are created by removing the gravitational effect of a tri-axial Moon model and the acceleration of the Sun and other planets from the raw Doppler tracking data. Cubic polynomials are fitted to the residuals for smoothing and estimation of accelerations. The accelerations are mapped to a constant orbit height of 100 km altitude above the Moon's surface. The point masses are introduced directly below the trajectory, with a depth of 50 km below the surface, and their magnitudes are estimated. Additional information is given in Wong et al. (1971), such as the restriction to 100 parameters in the estimation process due to implementation as well as the step-wise solutions in north–south bands, which usually cover 8 trajectories – with 48 elements in the estimated state vectors – and around 50 point masses below the tracks.

Figure 2 Modelling the gravity field of the asteroid 216 a.k.a. Kleopatra by point mass mascons according to Chanut et al. (2015). (a)  Polyhedron model of the asteroid; (b) geometrically sub-division of one tetrahedron by parallel planes with the mascon locations as large dots. Panel (b) was converted to grey-scale for this article.

2.3 Point mass mascons for irregular celestial bodies

Point mass mascons are also used in a different way to determine the gravity field of irregular celestial bodies. An example can be found in Chanut et al. (2015), where the gravity field of the asteroid 216 – also known as Kleopatra – is predicted by polyhedron models and point mass mascons. In the case of asteroids, the irregular shape is observed by optical instruments first; only in rare cases, orbiters investigate the gravity field directly. The observed shape is approximated by tetrahedrons with three corners on the surface and one in the geometrical centre of the asteroid (see Fig. 2a). Point masses are located then, either one per tetrahedron in its geometrical centre or three in the centres of a geometrically sub-divided tetrahedron (see Fig. 2b). Assuming a constant density of the asteroid and a known total mass, the mass per mascon is assigned to a value proportional to the surrounding volume, and the gravity field around the object can be predicted.

Properties

The point mass mascons have closed formulas for potential and explicit partial derivatives, which identify them as type A mascons in the threefold scheme.

+

The method is very easy to implement and requires only a few computational resources.

+

The gradient and all other field quantities are found without quadrature.

−

The model is singular for the potential and the gradient at the location of the point masses.

−

In the case of the lunar gravity field, assumptions are required for the location and depth below ground, as the Doppler tracking data and the observation geometry do not allow a detection of this information from the measurement.

It should be pointed out that the modelling by point masses is applied, for example, in Baur and Sneeuw (2011) or in Barthelmes (1986), Claessens et al. (2001), and Lin et al. (2014) without being labelled as a mascon approach by the authors, and that in the latter examples also the positions of the masses are estimated for regional studies of the Earth's gravity field. An iterative algorithm is developed and justified via quasi-orthogonality in the sense of an inner product in Barthelmes (1986). To stabilise the optimisation process, the possible movement per point mass shall be restricted in depth but also in radial or tangential direction with respect to an initial position.

2.4 Planar disc mascons

As a response to Muller and Sjogren (1968b), an article by Conel and Holstrom (1968) presented a physical interpretation of the ringed lunar maria, according to which former impact craters are filled afterwards by denser material. The authors experiment in the modelling of the mass anomalies with an arrangement of planar discs of finite thickness inside the impact craters and demonstrate a better post-fit to the residual Doppler tracking data for Mare Serenitatis.

The obvious issues of point masses are discussed in Wong et al. (1971):

• the singularities of the model at the centres;

• bad fitting of the residual tracking data in the equatorial zone of the Moon due to the observation geometry;

• and combination issues with the spherical harmonic models (of very low degree and order at the time).

To overcome these problems, finite mass elements are suggested for modelling the gravitational anomalies, which also agrees with the physical ideas in Conel and Holstrom (1968).

For a simple and efficient solution, the finite mass elements are chosen to be oblique rotational ellipsoids, also known as spheroids (Wong et al., 1971). The gravitational potential of a spheroid and its gradient are derived in Moulton (1960, pp. 119–132). On the one hand, the gravitational potential requires a series expression:

$$\begin{array}{}\text{(3)}& V={\displaystyle \frac{M}{R}}\left[\mathrm{1}+{\displaystyle \frac{b^{\mathrm{2}}}{\mathrm{10}}}{\displaystyle \frac{x_{\mathrm{P}}^{\mathrm{2}}+y_{\mathrm{P}}^{\mathrm{2}}-\mathrm{2}{z}_{\mathrm{P}}^{\mathrm{2}}}{R^{\mathrm{4}}}}{e}^{\mathrm{2}}+\mathcal{O}\left(e^{\mathrm{3}}\right)\right],\end{array}$$

with

• M: total mass of the spheroid (the gravitational constant is neglected in this exercise of the book);

• $R=\Vert {\mathit{r}}_{\mathrm{P}}\Vert$: Euclidean distance between the spheroid's centre and the calculation point; ${\mathit{r}}_{\mathrm{P}}=(x_{\mathrm{P}},y_{\mathrm{P}},z_{\mathrm{P}})$ outside the body,

• b: semi-minor axis of the spheroid (and semi-major axis a);

• $e=\sqrt{\frac{a^{\mathrm{2}}-b^{\mathrm{2}}}{a^{\mathrm{2}}}}$ numerical eccentricity.

On the other hand, the gradient of the potential can be derived in a closed formula. In Wong et al. (1971), the semi-minor axis b is then squeezed to zero, which leads to the attraction of a circular and planar disc. The article provides the gradient of a single disc in the form

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\ddot{x}& =-{\displaystyle \frac{\mathrm{3}Gm}{\mathrm{2}{a}^{\mathrm{3}}}}\left({\displaystyle \frac{-\sqrt{k}}{(\mathrm{1}+k)}}+\arcsin \left({\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+k}}}\right)\right)x\\ \ddot{y}& =-{\displaystyle \frac{\mathrm{3}Gm}{\mathrm{2}{a}^{\mathrm{3}}}}\left({\displaystyle \frac{-\sqrt{k}}{(\mathrm{1}+k)}}+\arcsin \left({\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+k}}}\right)\right)y\\ \ddot{z}& ={\displaystyle \frac{\mathrm{3}Gm}{a^{\mathrm{3}}}}\left({\displaystyle \frac{\mathrm{1}}{\sqrt{k}}}-\arcsin \left({\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}+k}}}\right)\right)z,\end{array}\end{array}$$

where k fulfils the quadratic equation

$$\begin{array}{}\text{(5)}& k^{\mathrm{2}}{a}^{\mathrm{2}}+(a^{\mathrm{2}}-\left(x^{\mathrm{2}}+y^{\mathrm{2}}+z^{\mathrm{2}}\right))k-z^{\mathrm{2}}=\mathrm{0}.\end{array}$$

To bring the expressions of the gradient in Moulton (1960) and Wong et al. (1971) into an analogous form, the identity $\arcsin \mathit{\zeta }=\arctan (\mathit{\zeta }/\sqrt{\mathrm{1}+{\mathit{\zeta }}^{\mathrm{2}}})$ must be kept in mind. It also turns out that the value k is linked to the numerical eccentricity of the spheroid by the relation $e=(\mathrm{1}/\sqrt{\mathrm{1}+k})$.

These planar disc mascons must be rotated and translated on the surface or close to it onto different locations, which is only implicitly indicated due to the definition of the coordinates $(x,y,z)$ with respect to the centre of each disc (Wong et al., 1971).

Properties

The planar disc mascons have closed formulas for explicit partial derivatives of the potential, which identifies them as type A mascons.

+

The closed formulas do not require any integration for the gradient.

+

The surface elements all have the same shape, size and area for each mascon.

−

The potential of a mascon requires a series expansion.

−

The model is singular for the potential at the centre of the disc.

−

The surface elements do not cover the complete surface, even in a global analysis.

−

Most points within a disc are either above or below the spherical surface.

In fact, the planar disc mascons are a kind of a simple layer potential, but without implicit or explicit integration for the gradient.

3 Simple layer potential and its regional subdivision

Modelling the gravitational potential by a simple layer was well known in geodesy and became popular around 1970.

The method can be applied to the complete potential or to a residual field after subtracting a reference field. The basic idea is to condensate the (remaining) in-homogeneous mass distribution onto the surface 𝒮, either the topography itself or a simpler reference like a sphere or spheroid (Koch and Witte, 1971; Morrison, 1971).

The gravitational potential of the layer is given by

$$\begin{array}{}\text{(6)}& V\left({\mathit{r}}_{\mathrm{P}}\right)=G\underset{\mathcal{S}}{\iint }{\displaystyle \frac{\mathit{\sigma }\left(\mathrm{\Omega }\right)}{\mathrm{\ell }\left(\mathrm{\Omega },{\mathit{r}}_{\mathrm{P}}\right)}}\mathrm{d}\mathrm{\Omega },\end{array}$$

with

• V(r_P): gravitational potential at the calculation point r_P;

• σ(Ω): location-dependent surface density;

• G: gravitational constant;

• $\mathrm{\ell }(\mathrm{\Omega },{\mathit{r}}_{\mathrm{P}})=\sqrt{(x_{\mathrm{P}}-x{)}^{\mathrm{2}}+(y_{\mathrm{P}}-y{)}^{\mathrm{2}}+(z_{\mathrm{P}}-z{)}^{\mathrm{2}}}$:
  Euclidean distance² between calculation point ${\mathit{r}}_{\mathrm{P}}=(x_{\mathrm{P}},y_{\mathrm{P}},z_{\mathrm{P}})$ and all surface points $\mathit{P}=(x,y,z)$, with P∈𝒮;

• dΩ: the differential surface element.

In the mascon version of the simple layer potential, the surface 𝒮 is sub-divided into smaller regions 𝒮_q – which are called surface elements or patches in this article – where the density is assumed to be constant. This leads to the mascon representation of the (residual) potential:

$$\begin{array}{}\text{(7)}& {\displaystyle }{V}_q\left({\mathit{r}}_{\mathrm{P}}\right)=G{\mathit{\sigma }}_q\underset{{\mathcal{S}}_q}{\iint }{\displaystyle \frac{\mathrm{1}}{\mathrm{\ell }(\mathrm{\Omega },{\mathit{r}}_{\mathrm{P}})}}\mathrm{d}\mathrm{\Omega },\text{(8)}& {\displaystyle }V\left({\mathit{r}}_{\mathrm{P}}\right)=\sum _{q=\mathrm{1}}^Q{V}_q\left({\mathit{r}}_{\mathrm{P}}\right).\end{array}$$

A linear combination Eq. (8) of all mascons – where the summation weights σ_q are included in the potential V_q(r_P) per mascon – generates the potential of the simple layer again. On the one hand, it should be pointed out that the method is applied, for example, in Koch and Witte (1971) without being labelled as a mascon approach. On the other hand, all the following mascon approaches are based on the simple layer potential with discrete surface elements and constant surface densities, which justifies the mascon label here.

3.1 Lumped spherical harmonics as mascons

Solving the Laplace equation in spherical coordinates leads to the spherical harmonic functions as a natural basis for gravity field modelling. An adequate linear combination of spherical harmonic functions can also be used to define localising base functions like the mascons in the spectral domain. Due to this combination over all degrees and orders, the result is sometimes labelled as the “lumped spherical harmonic approach” (Klosko et al., 2009).

Firstly, the gravity field is decomposed into a static field and its temporal variations:

$$\begin{array}{}\text{(9)}& V=V_{\mathrm{0}}+V_t.\end{array}$$

The static field and the mascons are represented by spherical harmonic synthesis. According to Heiskanen and Moritz (1967), Koch and Witte (1971) and Seeber (2003), the potential is given by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{V}_{\mathrm{0}}\left({\mathit{\lambda }}_{\mathrm{P}},{\mathit{\theta }}_{\mathrm{P}},r_{\mathrm{P}}\right)& ={\displaystyle \frac{GM}{r}}\sum _{l=\mathrm{0}}^L{\left({\displaystyle \frac{R}{r}}\right)}^l\sum _{m=\mathrm{0}}^l{\stackrel{\mathrm{‾}}{P}}_{l,m}\left(\cos {\mathit{\theta }}_{\mathrm{P}}\right)\\ & \times ({\stackrel{\mathrm{‾}}{C}}_{l,m}\cos m{\mathit{\lambda }}_{\mathrm{P}}+{\stackrel{\mathrm{‾}}{S}}_{l,m}\sin m{\mathit{\lambda }}_{\mathrm{P}}),\end{array}\end{array}$$

with

• $V_{\mathrm{0}}({\mathit{\lambda }}_{\mathrm{P}},{\mathit{\theta }}_{\mathrm{P}},r_{\mathrm{P}})$: potential of the static field;

• $({\mathit{\lambda }}_{\mathrm{P}},{\mathit{\theta }}_{\mathrm{P}},r_{\mathrm{P}})$: spherical coordinates of the evaluation point r_P, i.e. longitude λ_P, co-latitude θ_P and radius r_P;

• GM: product of gravitational constant G and the mass of the celestial body M;

• R: radius or semi-major axis of the spherical or ellipsoidal reference body;

• ${\stackrel{\mathrm{‾}}{P}}_{l,m}\left(\cos \mathit{\theta }\right)$: fully normalised Legendre functions;

• $\mathit{\{}{\stackrel{\mathrm{‾}}{C}}_{l,m},{\stackrel{\mathrm{‾}}{S}}_{l,m}\mathit{\}}$: fully normalised spherical harmonic coefficients, also known as Stokes coefficients.

The approach arose at GSFC when analysing the data of the GRACE mission, and it is presented in a sequence of articles (Rowlands et al., 2005; Lemoine et al., 2007; Klosko et al., 2009; Rowlands et al., 2010; Luthcke et al., 2013).

The mascons are generated in the spectral domain by (time-dependent) delta Stokes coefficients or differential Stokes coefficients of a simple layer:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathrm{\Delta }{\stackrel{\mathrm{‾}}{C}}_{l,m}^q\left(t\right)& ={\displaystyle \frac{\left(\mathrm{1}+k_{\mathrm{l}}^{\prime }\right)R^{\mathrm{2}}}{(\mathrm{2}l+\mathrm{1})M}}{\mathit{\sigma }}_q\left(t\right)\underset{{\mathcal{S}}_q}{\iint }{\stackrel{\mathrm{‾}}{P}}_{l,m}\left(\cos \mathit{\theta }\right)\cos m\mathit{\lambda }\mathrm{d}\mathrm{\Omega },\\ \mathrm{\Delta }{\stackrel{\mathrm{‾}}{S}}_{l,m}^q\left(t\right)& ={\displaystyle \frac{\left(\mathrm{1}+k_{\mathrm{l}}^{\prime }\right)R^{\mathrm{2}}}{(\mathrm{2}l+\mathrm{1})M}}{\mathit{\sigma }}_q\left(t\right)\underset{{\mathcal{S}}_q}{\iint }{\stackrel{\mathrm{‾}}{P}}_{l,m}\left(\cos \mathit{\theta }\right)\sin m\mathit{\lambda }\mathrm{d}\mathrm{\Omega },\end{array}\end{array}$$

with the Love numbers $k_{\mathrm{l}}^{\prime }$ for considering the loading effects of the extra masses on the surface.

The mascon solutions of the JPL are published online (https://earth.gsfc.nasa.gov/geo/data/grace-mascons, last access: 23 September 2022); the mascon solution of the CSR can also be found online (http://www.csr.utexas.edu/grace/, last access: 23 September 2022). As the formulas require standard techniques of geoscience, other groups are also working with these kind of mascons (e.g. Andrews et al., 2015; Krogh, 2011).

The lumped spherical harmonic approach can be used for any (almost spherical) body, but the approach is introduced for analysing the temporal variations of Earth's gravity field due to the variable water storage in particular. Taking into account that a uniform layer of 1 cm fresh water within an area of 1 m² has a mass of around 10 kg, the density is re-written in Rowlands et al. (2010) and Luthcke et al. (2013) as σ_q=10H_q – in Save et al. (2016) the factor σ_q=10.25H_q is used instead – to express the results in centimetres of equivalent water height. Each mascon is determined by a spherical harmonic synthesis

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{H}_q({\mathit{r}}_{\mathrm{P}},t)& ={\displaystyle \frac{M}{\mathrm{40}\mathit{\pi }{R}^{\mathrm{2}}}}\sum _{l=\mathrm{0}}^L\left({\displaystyle \frac{\mathrm{2}l+\mathrm{1}}{\mathrm{1}+k_{\mathrm{l}}^{\prime }}}\right)\sum _{m=\mathrm{0}}^l{\stackrel{\mathrm{‾}}{P}}_{l,m}\left(\cos {\mathit{\theta }}_{\mathrm{P}}\right)\\ & \times \left(\mathrm{\Delta }{\stackrel{\mathrm{‾}}{C}}_{l,m}^q\right(t)\cos m{\mathit{\lambda }}_{\mathrm{P}}+\mathrm{\Delta }{\stackrel{\mathrm{‾}}{S}}_{l,m}^q(t\left)\sin m{\mathit{\lambda }}_{\mathrm{P}}\right)\end{array}\end{array}$$

on the spherical surface r=R and with the upward continuation term $(R/r{)}^l$ in the synthesis formula, if necessary.

If the maximum degree L of the expansion is large enough, expression Eq. (12) forms a “two-dimensional step function” on the sphere 𝒮 (see Fig. 3), with

$$\begin{array}{}\text{(13)}& H_q({\mathit{r}}_{\mathrm{P}},t)=\left\{\begin{array}{ll}{H}_q& \text{in the region of interest, i. e. }{\mathcal{S}}_q\\ \mathrm{0}& \mathrm{outside}.\end{array}\right.\end{array}$$

A straightforward sub-division of a sphere is given by a longitude–latitude grid, i.e. all boundaries are either part of parallel circles or of meridians. In this case, the integrals Eq. (11) have the differential surface element dΩ = cos θdθdλ of the unit sphere, and the integration can be obtained by recursion formulas of integrated Legendre functions.

Figure 3 Sub-division in a longitude–latitude grid and a differential volume element of the sphere (Abedini et al., 2021a). In an alternative interpretation, the figure illustrates a two-dimensional step function on the sphere, which is either zero (white) outside and constant but non-zero (grey) inside a mascon element.

The size and shape of the surface elements vary between the publications:

• Lemoine et al. (2007) and Rowlands et al. (2005) present a separation of the region of interest into surface elements of equal angles with the dimension $\mathrm{4}{}^{\circ }\times \mathrm{4}{}^{\circ }$, while Krogh (2011) defines patches of the dimensions $\mathrm{1.25}{}^{\circ }\times \mathrm{1.5}{}^{\circ }$ and $\mathrm{1.5}{}^{\circ }\times \mathrm{1.5}{}^{\circ }$.

• Equal areas within a longitude–latitude grid can be obtained by stretching or shrinking one of the angles dependent on latitude, which is discussed already in Morrison (1971) and applied in experiments of Rowlands et al. (2010) and Andrews et al. (2015).

• In Klosko et al. (2009), the surface elements have – at least in the corresponding Fig. 4 – more complex boundaries. The lines are still along parallel circles and meridians, but they are combined in such a way that the mascon patches fill irregular shapes of sub-basins within the Mississippi basin.

• In the CSR solution, the equal area per mascon is considered to be more relevant than a simple sub-division or a complete coverage of the sphere (Save et al., 2016). A geodesic grid with 40 962 vertices is generated by iteration, and the mascon patches are located in the centres. The shapes of the patches are either hexagonal or pentagonal, and the elements cover approximately equal areas of around 1^∘ diameter.

The regularisation techniques for equiangular patches are discussed in Abedini et al. (2021b) for other types of mascons, but the final recommendation to consider herein the area size should be transferable to the lumped spherical harmonic approach as well.

Figure 4 Definition of mascon surface elements in the Mississippi basin. The figure originates from Klosko et al. (2009), but it is converted to grey-scale for this article.

Observation of GRACE

The mascons were introduced for analysing the Earth's gravity field in the mission GRACE (Gravity Recovery And Climate Experiment). The mission consisted of two identical satellites, which were launched in 2002 in a cooperation between NASA/JPL and the German DLR. The satellites fell around the Earth in one common and almost circular orbit, with a low altitude of originally 500 km in height. The positions were quasi-permanently observed by GPS receivers with three antennas, and onboard accelerometers with three axes were measuring the combined influence of all non-gravitational effects. The main observable was the variation of the distance between the two GRACE satellites, measured by microwaves in the K-band and Ka-band via a range-rate measurement system. The distance of ρ≈250 km between the satellite centres varied due to mass variations below, and the K-band provided a nominal accuracy of 10 µm for the range ρ and 0.5 µm s⁻¹ for the range-rate $\dot{\mathit{\rho }}$ (Seeber, 2003; Tapley et al., 2004).

The orbit observations and the gravity field parameters can be linked in different ways – e.g. the variational equation, the energy balance approach, the short arc approach or the acceleration approach – which are sketched, for example, in Liu (2008). The details are not in the focus of this work, but most methods require the gradient of the gravitational potential again.

Gradient of the lumped spherical harmonic mascons

For the range-rate $\dot{\mathit{\rho }}$ in the lumped harmonic approach, the relationship is found in Luthcke et al. (2013) by the chain rule

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial \dot{\mathit{\rho }}}{\partial H_q}}=\sum _{l=\mathrm{0}}^L\sum _{m=\mathrm{0}}^l{\displaystyle \frac{\partial \dot{\mathit{\rho }}}{\partial {\stackrel{\mathrm{‾}}{C}}_{l,m}}}{\displaystyle \frac{\partial \mathrm{\Delta }{\stackrel{\mathrm{‾}}{C}}_{l,m}^q}{\partial H_q}}+{\displaystyle \frac{\partial \dot{\mathit{\rho }}}{\partial {\stackrel{\mathrm{‾}}{S}}_{l,m}}}{\displaystyle \frac{\partial \mathrm{\Delta }{\stackrel{\mathrm{‾}}{S}}_{l,m}^q}{\partial H_q}}\end{array}$$

and analogously for range-acceleration $\ddot{\mathit{\rho }}$. The derivatives $\left\{\frac{\partial \mathrm{\Delta }{\stackrel{\mathrm{‾}}{C}}_{l,m}^q}{\partial H_q},\frac{\partial \mathrm{\Delta }{\stackrel{\mathrm{‾}}{S}}_{l,m}^q}{\partial H_q}\right\}$ are straightforward, as the formulas Eq. (11) are linear with respect to the surface density σ_q or the water height H_q.

Properties

The lumped spherical harmonic approach is a representative of the type B mascons.

+

The method is very easy to implement after a previous analysis of the GRACE observations by spherical harmonic functions.

+

After the determination of all delta Stokes coefficients, all other field quantities can be calculated by standard methods of spherical harmonic synthesis.

+

The required integration Eq. (11) can be solved by well-known recursion formulas or by numerical quadrature.

−

A high degree L of expansion might be required for straight boundaries and constant values within the two-dimensional step functions.

3.2 Spherical cap mascons

The planar disc mascon approach (in Sect. 2.4) is not satisfying from a geometrical view point, as most points within the element are either above or below the spherical surface. This can be avoided by introducing spherical caps instead of planar discs. Monthly solutions in terms of spherical cap mascons are calculated at the JPL, and the details can be found in Watkins et al. (2015).

To reduce the effort of quadrature for the gradient expression, a local mascon coordinate system is introduced for each element by rotation, where the centre of the mascon is equal to the new North Pole of the system. The new coordinates are the spherical distance γ and the azimuth ξ in the calculation point. The potential is still based on the simple layer theory leading to the integral

$$\begin{array}{}\text{(15)}& {\stackrel{\mathrm{‾}}{V}}_q\left({\mathit{r}}_{\mathrm{P}}\right)=R^{\mathrm{2}}{\mathit{\sigma }}_q\underset{\mathrm{0}}{\overset{\mathit{\alpha }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}{\displaystyle \frac{\mathrm{d}\mathit{\xi }\sin \mathit{\gamma }\mathrm{d}\mathit{\gamma }}{\mathrm{\ell }\left(\mathrm{\Omega },{\mathit{r}}_{\mathrm{P}}\right)}}\end{array}$$

for the potential of a spherical cap, with

• ${\stackrel{\mathrm{‾}}{V}}_q\left({\mathit{r}}_{\mathrm{P}}\right)$: gravitational potential per mascon in the local mascon coordinate system (the over-bar is introduced here to emphasise the rotated coordinate system);

• σ_q: product of the gravitational constant G and mass per mascon m_q divided by the area of the spherical cap, i.e. ${\mathit{\sigma }}_q=\frac{G{m}_q}{\mathrm{2}\mathit{\pi }(\mathrm{1}-\cos \mathit{\alpha })R^{\mathrm{2}}}$;

• R: radius of the spherical reference model;

• ℓ(Ω,r_P): distance between calculation point ${\mathit{r}}_{\mathrm{P}}=(x_{\mathrm{P}},y_{\mathrm{P}},z_{\mathrm{P}})$ and the points $\mathit{P}=(x,y,z)$ in the spherical cap (in the original paper, the Euclidean distance is noted down by d);

• α: radius of the spherical cap in radians.

Gradient of the spherical cap mascons

The gradient of the potential ${\stackrel{\mathrm{‾}}{V}}_q\left(r_{\mathrm{P}}\right)$ is calculated per mascon and rotated at the end to the original coordinate system. The iterated integration over the spherical distance and azimuth is reduced to a single integral by expressing the azimuthal component via elliptic integrals.

In the local mascon coordinate system, the gradient operator is of the form

$$\begin{array}{}\text{(16)}& \mathrm{\nabla }{\stackrel{\mathrm{‾}}{V}}_q={\left(\begin{array}{ccc}\frac{\partial {\stackrel{\mathrm{‾}}{V}}_q}{\partial r},& \mathrm{0},& \frac{\mathrm{1}}{r}\frac{\partial {\stackrel{\mathrm{‾}}{V}}_q}{\partial \mathit{\theta }}\end{array}\right)}^{\top },\end{array}$$

where θ is the spherical distance between the calculation point and the mascon centre. The formulas of the gradient of a spherical cap are derived in an inter-office memorandum at the JPL (R. Sunseri, unpublished data: Mass concentration modelled as a spherical cap 343R-11-00) – which is not available to us – and the results are quoted by Watkins et al. (2015):

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\displaystyle \frac{\partial V_q}{\partial r}}& =-{\mathit{\sigma }}_q{t}^{\mathrm{3}}\left({\displaystyle \frac{I_{\mathrm{2}}}{t}}-\cos \mathit{\theta }\cdot I_{\mathrm{1}}-\sin \mathit{\theta }\cdot I_{\mathrm{3}}\right)\\ {\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial V_q}{\partial \mathit{\theta }}}& =-{\mathit{\sigma }}_q{t}^{\mathrm{3}}\left(\sin \mathit{\theta }\cdot I_{\mathrm{1}}-\cos \mathit{\theta }\cdot I_{\mathrm{3}}\right),\end{array}\end{array}$$

with the abbreviation $t=\frac{R}{r}$ and the three integrals $\mathit{\{}{I}_{\mathrm{1}},I_{\mathrm{2}},I_{\mathrm{3}}\mathit{\}}$. The solution of the later ones requires complete elliptic integrals – first kind E(k) and second kind K(k) – and numerical integration in the spherical distance direction:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{I}_{\mathrm{1}}& =\int \sin \mathit{\gamma }\cos \mathit{\gamma }\left[{\displaystyle \frac{m^{\prime }}{\sqrt{l^{\prime }+\mathrm{1}}(l^{\prime }+\mathrm{1})}}E\left(k\right)\right]\mathrm{d}\mathit{\gamma }\\ I_{\mathrm{2}}& =\int \sin \mathit{\gamma }\left[{\displaystyle \frac{m^{\prime }}{\sqrt{l^{\prime }+\mathrm{1}}(l^{\prime }+\mathrm{1})}}E\left(k\right)\right]\mathrm{d}\mathit{\gamma }\\ I_{\mathrm{3}}& =\int \sin \mathit{\gamma }\cos \mathit{\gamma }\left[{\displaystyle \frac{m^{\prime }\left(E\right(k)-(\mathrm{1}-l^{\prime }\left)K\right(k\left)\right)}{\sqrt{l^{\prime }+\mathrm{1}}(l^{\prime }+\mathrm{1})l^{\prime }}}\right]\mathrm{d}\mathit{\gamma },\end{array}\end{array}$$

with the auxiliary expressions

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& n=\mathrm{1}+t^{\mathrm{2}}-\mathrm{2}t\cos \mathit{\theta }\cos \mathit{\gamma }\\ & m^{\prime }=\mathrm{4}/n^{\mathrm{3}/\mathrm{2}}\\ & l^{\prime }=\mathrm{2}t\sin \mathit{\theta }\sin \mathit{\gamma }/n\\ & k^{\mathrm{2}}=\mathrm{2}{l}^{\prime }/(l^{\prime }+\mathrm{1}).\end{array}\end{array}$$

Properties

The mascon potential is calculated by quadrature, and analytical derivatives have been derived, which leads to a class A mascon in the threefold scheme.

+

The two-dimensional quadrature for the gradients are reduced to a one-dimensional integration.

+

The calculation takes place only in the spatial domain and avoids the truncation error of spherical harmonic synthesis.

+

The surface elements have all the same shape, size and area for each mascon.

−

The surface elements do not cover the complete surface, even in a global analysis.

−

The model is singular for the potential and the gradient at the location of the centre of the spherical cap.

−

A straightforward implementation of the formulas Eq. (19) also leads to an undefined expression when the calculation point is identical to the centre of the spherical cap. One finds then that t=1 and θ=0, and in consequence $l^{\prime }=\mathrm{0}/\mathrm{0}$. A solution might be given in the unavailable inter-office memorandum.

3.3 Mascons via quadrature of the simple layer potential

To avoid truncation errors and aliasing into coefficients of lower degree and order via the spherical harmonic expansion Eq. (14), a complete numerical integration is suggested in Abedini et al. (2021a, b). The potential is represented – in our notation of Sect. 3 – by the formula

$$\begin{array}{}\text{(20)}& T\left({\mathit{r}}_{\mathrm{P}}\right)=-G\underset{\mathcal{S}}{\iint }{\displaystyle \frac{\mathit{\sigma }\left(\mathrm{\Omega }\right)}{\mathrm{\ell }\left(\mathrm{\Omega },{\mathit{r}}_{\mathrm{P}}\right)}}\mathrm{d}\mathrm{\Omega }\end{array}$$

and is evaluated by numerical quadrature when necessary. The extra minus sign was likely introduced by the authors due to non-geodetic literature, as physical textbooks often use the definition $\ddot{\mathit{r}}=-\mathrm{\nabla }V$.

The derivatives of the range-rate $\dot{\mathit{\rho }}$ with respect to the surface density are quasi-decomposed by the chain rule

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial \dot{\mathit{\rho }}}{\partial {\mathit{\sigma }}_q}}={\displaystyle \frac{\partial \dot{\mathit{\rho }}}{\partial \mathit{X}}}{\displaystyle \frac{\partial \mathit{X}}{\partial {\mathit{\sigma }}_q}}+{\displaystyle \frac{\partial \dot{\mathit{\rho }}}{\partial \dot{\mathit{X}}}}{\displaystyle \frac{\partial \dot{\mathit{X}}}{\partial {\mathit{\sigma }}_q}}\end{array}$$

into geometrical components $\left\{\frac{\partial \dot{\mathit{\rho }}}{\partial \mathit{X}},\frac{\partial \dot{\mathit{\rho }}}{\partial \dot{\mathit{X}}}\right\}$ and dynamical components $\left\{\frac{\partial \mathit{X}}{\partial {\mathit{\sigma }}_q},\frac{\partial \dot{\mathit{X}}}{\partial {\mathit{\sigma }}_q}\right\}$, with

• $\mathit{X}={\mathit{X}}_{\mathrm{2}}-{\mathit{X}}_{\mathrm{1}}$: difference vector between the satellites' centre positions;

• $\dot{\mathit{X}}={\dot{\mathit{X}}}_{\mathrm{2}}-{\dot{\mathit{X}}}_{\mathrm{1}}$: difference vector between the satellites' centre velocities.

As the range-rate $\dot{\mathit{\rho }}$ is calculated by $\dot{\mathit{\rho }}=\frac{{\mathit{X}}^{\top }\dot{\mathit{X}}}{\Vert \mathit{X}\Vert }$, the geometrical components are known and can be differentiated with respect to positions and velocities.

The dynamical components are determined by the variational equation

$$\begin{array}{}\text{(22)}& \ddot{\mathit{\xi }}={\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial t^{\mathrm{2}}}}\left\{{\displaystyle \frac{\partial \mathit{X}}{\partial {\mathit{\sigma }}_q}}\right\}={\mathrm{\nabla }}^{\mathrm{2}}U\left(\mathit{X}\right){\displaystyle \frac{\partial \mathit{X}}{\partial {\mathit{\sigma }}_q}}+{\displaystyle \frac{\partial \mathrm{\nabla }T}{\partial {\mathit{\sigma }}_q}},\end{array}$$

with $U=\frac{GM}{\Vert {\mathit{r}}_{\mathrm{P}}\Vert }$ being the potential of the Kepler problem. The equation is solved similarly to an orbit integration, with the initial values ξ=0 and $\dot{\mathit{\xi }}=\mathbf{0}$ for each arc, each satellite and each mascon. The integration error is limited by applying the method only to short arcs over the region of interest, e.g. Greenland in Abedini et al. (2021a).

Properties

The approach does not fit into the threefold scheme.

+

The method avoids truncation errors and aliasing by integration in the spatial domain.

+

The surface elements cover the complete surface in a global analysis.

−

The potential and the gradient require numerical quadrature.

−

The variational equations lead to a high computational burden, which is already admitted in Abedini et al. (2021a).

4 Mascons as a post-processing tool

Since the successful GRACE mission, it is also possible to observe the temporal variations of the gravity field. The standard output of these investigations are monthly solutions of spherical harmonic coefficients, which are meanwhile complemented by mascon solutions in the same time span by several research centres.

The question arises whether it is possible to estimate local variations from the spherical harmonic solutions by post processing. This is of particular interest for the ice masses and glaciers in Greenland, Antarctica and Alaska as well as for the highly variable water masses in the large water basins, which dominate the time-variable part of the gravity field.

The spherical harmonic functions have a global support, which contradicts a regional analysis. Another problem is the noise in the coefficients, which is overcome by filtering and de-striping techniques at the cost of the spatial resolution. To estimate regional mass changes, it can be helpful to determine an adequate field quantity by spherical harmonic synthesis and to analyse this newly generated signal by another base function with local support (Ran et al., 2018).

4.1 Spherical cap mascon as a post-processing tool

Schrama et al. (2014) use the term mascon for post processing of a time series of Stokes coefficients $\left\{{\stackrel{\mathrm{‾}}{C}}_{l,m}\left(t\right),{\stackrel{\mathrm{‾}}{S}}_{l,m}\left(t\right)\right\}$. The goal is the determination of local mass variations in the ice shields and glaciers based on a time series of spherical harmonic coefficients.

A long-term mean value $\left\{\langle {\stackrel{\mathrm{‾}}{C}}_{l,m}\rangle ,\langle {\stackrel{\mathrm{‾}}{S}}_{l,m}\rangle \right\}$ per coefficient is calculated and subtracted, and a Gauß-filter ${W_{\mathrm{l}}}^{\mathrm{G}}$ in the spectral domain is applied by multiplication with the Stokes coefficients. To represent the equivalent water height instead of the potential, further standard factors are introduced:

$$\begin{array}{}\text{(23)}& \left\{\begin{array}{c}{c}_{l,m}^{\mathrm{w}}\left(t\right)\\ s_{l,m}^{\mathrm{w}}\left(t\right)\end{array}\right\}={\displaystyle \frac{a_{\mathrm{e}}{\mathit{\rho }}_{\mathrm{e}}(\mathrm{2}l+\mathrm{1})}{\mathrm{3}{\mathit{\rho }}_{\mathrm{w}}\left(\mathrm{1}+k_{\mathrm{l}}^{\prime }\right)}}{W_{\mathrm{l}}}^{\mathrm{G}}\left\{\begin{array}{c}{\stackrel{\mathrm{‾}}{C}}_{l,m}\left(t\right)-\langle {\stackrel{\mathrm{‾}}{C}}_{l,m}\rangle \\ {\stackrel{\mathrm{‾}}{S}}_{l,m}\left(t\right)-\langle {\stackrel{\mathrm{‾}}{S}}_{l,m}\rangle \end{array}\right\},\end{array}$$

with

• a_e: equatorial radius of the ellipsoidal Earth;

• ρ_e: mean density of the Earth;

• ρ_w: density of water;

• $k_{\mathrm{l}}^{\prime }$: Love numbers.

The spherical harmonic synthesis

$$\begin{array}{}\text{(24)}& \begin{array}{rl}h\left({\mathit{\lambda }}_{\mathrm{P}},{\mathit{\theta }}_{\mathrm{P}},t\right)& =\sum _{l=\mathrm{0}}^L\sum _{m=\mathrm{0}}^l{\stackrel{\mathrm{‾}}{P}}_{l,m}\left(\cos {\mathit{\theta }}_{\mathrm{P}}\right)\\ & \times \left(c_{l,m}^{\mathrm{w}}\right(t)\cos m{\mathit{\lambda }}_{\mathrm{P}}+s_{l,m}^{\mathrm{w}}(t\left)\sin m{\mathit{\lambda }}_{\mathrm{P}}\right)\end{array}\end{array}$$

provides the mass variations with respect to a long-term mean on a spherical surface. The equivalent water height $h({\mathit{\lambda }}_{\mathrm{P}},{\mathit{\theta }}_{\mathrm{P}},t)$ is then analysed by a set of localising base functions

$$\begin{array}{}\text{(25)}& h\left({\mathit{\lambda }}_{\mathrm{P}},{\mathit{\theta }}_{\mathrm{P}},t\right)=\sum _{q=\mathrm{1}}^Q{\mathit{\alpha }}_q\left(t\right){\mathit{\beta }}_q\left({\mathit{\psi }}_q,L,R\right)\end{array}$$

via least-squares estimation and the determination of the weights α_q(t). Each base function ${\mathit{\beta }}_q({\mathit{\psi }}_k,L,R)$ has the form

$$\begin{array}{}\text{(26)}& {\displaystyle }{\mathit{\beta }}_q=\sum _{l=\mathrm{0}}^L{\mathit{\gamma }}_{\mathrm{l}}\left(R\right){\stackrel{\mathrm{‾}}{P}}_{\mathrm{l}}\left(\cos {\mathit{\psi }}_q\right)\text{(27)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{l}}\left(R\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{R}{\int }}{\stackrel{\mathrm{‾}}{P}}_{\mathrm{l}}\left(\cos \mathit{\mu }\right)\sin \mathit{\mu }\mathrm{d}\mathit{\mu },\end{array}$$

which is equivalent to a spherical cap with the radius R, the maximum expansion degree L in the spectral domain and the location (λ_q,θ_q) for its centre.

Properties

Original GRACE data are not required, as the method is applied to the previous solution by spherical harmonics, which leads to type C mascons.

+

The estimation of the weights is a straightforward process via least-squares estimation.

+

The surface elements all have the same shape, size and area.

+

The required integration Eq. (27) can be solved by recursion formulas or numerical quadrature.

−

The effect of Gauß filtering and the temporal average on the solution's quality are difficult to predict; also, the chosen sampling in the spherical harmonic synthesis might have an effect on the estimated masses.

−

The surface elements do not cover the complete surface, even in a global analysis.

4.2 Point mass mascons as a post-processing tool

Ran (2017) and Ran et al. (2018) extend an idea of Baur and Sneeuw (2011) by combining point masses and the simple layer potential. The goal of the work is the estimation of mass variations over Greenland based on spherical harmonic coefficients. The GRACE solutions are used to derive the radial component of the gradient but with loading compensation in orbit altitude:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}\mathit{\delta }g\left({\mathit{r}}_{\mathrm{P}}\right)& =-{\displaystyle \frac{\partial V}{\partial r}}\\ & =-{\displaystyle \frac{GM}{r_{\mathrm{P}}^{\mathrm{2}}}}\sum _{l=\mathrm{1}}^L{\displaystyle \frac{l+\mathrm{1}}{\mathrm{1}+k_{\mathrm{l}}^{\prime }}}{\left({\displaystyle \frac{a}{r_{\mathrm{P}}}}\right)}^l\sum _{m=\mathrm{0}}^l{\stackrel{\mathrm{‾}}{P}}_{l,m}\left(\cos \mathit{\theta }\right)\\ & \times (\mathrm{\Delta }{\stackrel{\mathrm{‾}}{C}}_{l,m}\cos m{\mathit{\lambda }}_{\mathrm{P}}+\mathrm{\Delta }{\stackrel{\mathrm{‾}}{S}}_{l,m}\sin m{\mathit{\lambda }}_{\mathrm{P}}).\end{array}\end{array}$$

This signal is analysed – by least squares estimation of the surface densities ρ_q – by the simple layer in the region of interest:

$$\begin{array}{}\text{(29)}& \begin{array}{rl}\mathit{\delta }{q}_{\mathrm{P}}& =-{\displaystyle \frac{\partial }{\partial r}}\left\{G\sum _{q=\mathrm{1}}^Q{\mathit{\rho }}_q\iint {\displaystyle \frac{\mathrm{d}s}{\mathrm{\ell }(\mathrm{\Omega },r_{\mathrm{P}})}}\right\}=\\ & =-{\displaystyle \frac{\partial }{\partial r}}\left\{G\sum _{q=\mathrm{1}}^Q{\mathit{\rho }}_q{I}_{q,p}\right\}.\end{array}\end{array}$$

The integral I_q,p is approximated by quadrature, which evaluates the distances only in the nodes of a grid:

$$\begin{array}{}\text{(30)}& I_{q,p}=\iint {\displaystyle \frac{\mathrm{d}s}{\mathrm{\ell }\left(\mathrm{\Omega },{\mathit{r}}_{\mathrm{P}}\right)}}\approx \sum _{j=\mathrm{1}}^{K_q}{w}_{q,j}{\displaystyle \frac{\mathrm{1}}{l_{q,j,p}}},\end{array}$$

with

• $w_{q,j}=S_i/K_q$: weighting of the evaluated points;

• $l_{q,j,p}$: distance between the nodes and the evaluation point;

• S_q: the surface area of the mascon with the index q.

The Euclidean distance is expressed in spherical coordinates

$$l_{q,j,p}=\sqrt{r_{q,j}^{\mathrm{2}}+r_p^{\mathrm{2}}-\mathrm{2}{r}_{q,j}{r}_p\cos {\mathrm{\Psi }}_{q,j,p}},$$

with

• $r_p=\Vert {\mathit{r}}_{\mathrm{P}}\Vert$: distance of calculation point to the origin;

• $r_{q,j}=\Vert {\mathit{r}}_{q,j}\Vert$: distance of nodes within the patch to the origin;

• ${\mathrm{\Psi }}_{q,j,p}$: angle between the vectors r_q,j and r_P.

The observable of the study is then given by

$$\begin{array}{}\text{(31)}& \mathit{\delta }{q}_{\mathrm{P}}\approx G\sum _{q=\mathrm{1}}^Q{\mathit{\rho }}_q\sum _{j=\mathrm{1}}^{K_q}{w}_{q,j}{\displaystyle \frac{\left(r_{q,j}-r_p\cos {\mathrm{\Psi }}_{q,j,p}\right)}{{l_{q,j,p}}^{\mathrm{3}}}}.\end{array}$$

Properties

The method is applied to the previous solution by spherical harmonics, which leads to type C mascons.

+

The method is very easy to implement.

+

Integration per mascon element is replaced by a weighted sum of point masses located on a grid.

−

The model is singular for potential and radial derivatives at the location of the nodes.

−

Finding the weighting $w_{q,j}=S_i/K_q$ might be challenging for irregularly shaped patches.

5 Summary – mascons in gravity field modelling

Point mass models are an important tool for gravity field modelling due to their simplicity and efficiency. The point mass representation is used for celestial bodies with irregular shapes but also for Earth or the Moon on regional and global scales. Point mass mascons are also a key aspect of converting spherical harmonic solutions into regional mass variations, which supports the interpretation of geophysical processes.

Mascons represented by finite surface elements are based on the simple layer potential. These models form a subset of localising base functions for gravity field modelling. Without neighbourhood conditions, a solution close to the ground generates a discontinuous field. The discontinuity problem is damped for higher evaluation altitude or small patches. The constant density per mascon simplifies the interpretation of mass variations in comparison to other localising base functions (e.g. wavelets or Slepian functions), which vary within their region of interest. Planar disc and spherical cap mascons are radial symmetric base functions, while the other mascon concepts allow for patches with arbitrary shapes. In particular, the shape can consider the geometry of water basins, which reduces the leakage of signals in hydro-geodesy.