Gravitational inversion for arbitrary surface mass distribution—application to earth's time-variable gravity with topography

SUMMARY

1. INTRODUCTION

Since its launch in 2002, the dual-satellite mission of GRACE (Gravity Recovery And Climate Experiment) and its Follow-on have revolutionized the way we observe and understand the Earth's mass transport processes in the era of global changes (Tapley et al. 2004; see references in, e.g. https://grace.jpl.nasa.gov/publications/). The mass transports occur in Earth's surface fluids of atmosphere, oceans, cryosphere and land hydrology, as well as in solid Earth's interior of the crust, mantle and cores.

Over the years a progression of the theoretical formulation of planetary gravity, or more specifically the time-variable gravity (TVG) produced by an arbitrary mass redistribution, has been developed and utilized in pace with the continual improvements and higher observation and modelling demands for the GRACE research. The 3-D general formula was developed in fact far prior to GRACE (e.g. Chao & Gross 1987; Chao et al. 1987). Its 2-D reduced version soon took over in practice not surprisingly because most targeted mass transports do occur on Earth's surface, and the spherical-Earth approximation sufficed in all early applications. For GRACE, that progressively led to the convenient outcome in the form of equivalent water thickness (EWT) or surface ‘mascon’.

This course of progression has unintentionally fallen in a contorted path that raises probable misconceptions. In particular, the spherical-surface special case was often misunderstood in the literature to be taken as the ‘standard’ form (e.g. when referencing to Wahr et al. 1998), on top of which all applications were to stem. Furthermore, later pursuits for higher precision were mis-considered as ‘corrections’ and often treated as extensions from the ‘standard’ or even as something extraneous. Meanwhile the associated issue of non-uniqueness in the gravitational inversion, one that can be traced back to the work of Stokes (1867), came to bear but drew inadequate consideration even in cases that ought to.

Increasing demands in recent years for more precise modelling called into attention the non-spherical shape of the Earth, whose gravitational effects surely do get felt by the orbiting GRACE satellites. Recent work based on various approaches of approximation (e.g. Li et al. 2017; Ditmar 2018; Ghobadi-Far et al. 2019; Yang et al. 2022) aimed for refined GRACE mascon solutions that can account for the Earth surface shape, including the ellipsoidal figure and the topography, whereof the actual surface mass redistributions reside.

In this paper we utilize the multipole formalism (e.g. Chao & Shih 2021) to reinstate the logic and give the exact and general equations for TVG of the planet Earth. The formulation proves the uniqueness of the gravitational inversion of an arbitrary 2-D surface, which leads to a mathematically rigorous algorithm for determining the surface mass distribution in terms of spherical harmonics (SH) up to a truncated finite degree.

2. GRAVITATION OF A GENERAL 3-D BODY

2.1 Multipole expansion

We first summarize the multipole formalism for the Newtonian gravitation, the details of which was given in Chao & Gross (1987); a more general review was given by Chao & Shih (2021). The gravitational field U produced externally by a 3-D central body (a ‘potato’) with mass density distribution ρ(r₀) (suffix 0 referring to the body) is:

where G is the gravitational constant. The addition theorem |$1/| {{{\bf r}}{{{\bf r}}}_0} |{\rm{\ }} = \mathop \sum_{n = 0}^\infty \mathop \sum _{m = - n}^n \frac{{4\pi \ {r}_0^n}}{{( {2n + 1} )\ {r}^{n + 1}}}{Y}_{nm}( {{\Omega }_0} )Y_{nm}^*( \Omega )\ $| then leads to the multipole expansion:

where r = (radius r, co-latitude θ, East longitude λ) locates the field point, Ω is the solid angle representing θ and λ so that dV = r²dr dΩ and dΩ = sinθ dθ dλ and * denotes complex conjugation. The fully normalized complex surface SH function |${Y}_{nm}( \Omega )\ = \ ( { - 1{)}^m{{( {\frac{{( {2n + 1} )( {n - m} )!}}{{4\pi ( {n + m} )!}}} )}}^{1/2}{P}_{nm}( {\cos \theta } )\exp (im\lambda } )$| (P_nm denoting the associated Legendre function) of degree n ( = 0, 1,…, ∞) and order m ( = −n,..0,…,n), henceforth abbreviated as [n, m], is subject to the orthonormality:

For convenience, we stipulate that our geocentric coordinate origin coincides with the centre of mass of the body. The nominal co-latitude θ and longitude λ refer to mathematical coordinates not necessarily having any terrestrial connotation, so that the practical issue of mapping the artificial ‘grid’ points among different types of Earth surface models (such as spheroidal) need not enter.

The complex-valued volume integral in the bracket of eqs (2), namely (dropping the suffix 0):

is the [n, m]th SH multipole moments of the exterior type belonging to the body of mass density ρ(r) (Chao & Shih 2021), henceforth simply the [n, m]th multipoles. The n = 0 monopole is Q₀₀ = (4π)^−1/2*(total mass); the three n = 1 dipoles signify the location of the centre of mass so they vanish under the above stipulation about the coordinate origin. The gravitational ‘anomalies’ start with n = 2 quadrupoles. Expression (2) is real-valued, as it should, by virtue of the symmetry Y_n-m = (−1)^mY_nm*. A counterpart set of formulae in terms of multipoles of the interior type applies to the interior gravitational field (for details see Chao & Shih 2021); they are not relevant here.

Now compare eq. (4) with the conventional expansion of U in real-valued form in planetary geodesy (e.g. Kaula 1966; Torge 1989) |$U( {{\bf r}} )\ = \frac{{GM}}{a}\ \mathop \sum _{n\ = \ 0}^\infty \mathop \sum_{m\ = \ 0}^n {( {\frac{a}{r}} )}^{n + 1}{\tilde{P}}_{nm}( {\cos \theta } )( {{C}_{nm}\cos m\lambda + {S}_{nm}\sin m\lambda } )$|⁠, where M is the total mass of the planet and a its equatorial radius, the 4π-normalized |${\tilde{P}}_{nm}$|= [(2−δ_m0)(2n + 1)(n − m)!/(n + m)!]^1/2P_nm, and C_nm and S_nm are the [n, m]th dimensionless Stokes coefficients which comprise the GRACE Level-2 data. We see readily that the Stokes coefficients are just normalized multipoles (Chao & Gross 1987):

Note that the presence of a here is just for scaling, not implying any connection with certain sphere per se.

2.2 Time-variable gravity

Gravity is generated dependent on the ‘state’ mass distribution at any given moment, whether or not the mass distribution is static or time-variable. Here we adopt the Eulerian (as opposed to Lagrangian) description for the TVG field in the form of time-variable Δρ residing at a given location in space. Definition (4) leads immediately to:

for a time-variable mass redistribution. The Δ enters likewise in relation (5) on both sides.

2.3 Non-uniqueness in 3-D gravitational inversion

The above constitutes the forward problem, where the gravity (change) is calculable directly from the source mass (re)distribution. The gravitational inverse problem by contrast refers to the reconstruction of the source mass density (re)distribution basing on the known gravity (change). The well-known non-uniqueness of the 3-D gravitational inversion (e.g. Stokes 1867; Parker 1994) has been elucidated from the perspective of multipoles by Chao (2005). (Incidentally, the present uniqueness problem is not to be confused with the uniqueness theorem for a Laplace field given its boundary conditions.) Specifically, there are infinitely many possible mass distribution that have the same set of multipoles. Knowing all the multipoles or Stokes coefficients of the body, even if precisely and completely out to degree infinity, is grossly insufficient to uniquely determine the mass distribution that generates the gravitational field in the first place.

Mathematically, the generic gravitational inversion problem can be cast in the framework of the linear Hilbert function space with the notion of inner product. We do so in Appendix A, recognizing that a multipole (4) is the inner product (or ‘projection’) of the density distribution ρ(r) onto the regular solid harmonic rⁿY_nm(Ω), or in shorthand 〈ρ, rⁿY_nm〉 (see Appendix A) using the bra-ket notation that prevails (in fact is indispensable) in quantum mechanics (e.g. Dirac 1958). The root of the present non-uniqueness lies in the fact that the basis functions {rⁿY_nm}, obeying Laplace's equation, fail to form a complete set to span the Hilbert space H in 3-D; they are complete only for the sub-space of H consisting of Laplace fields, not other functions in general (case in point ρ(r)). Consequently the corresponding set of projections of ρ(r) onto that incomplete set of basis functions, namely the multipoles, provides only partial information insufficient to allow a full reconstruction of ρ(r) itself; the ‘degree of deficiency’ is as much as n(n−1)/2 for each degree n (Chao 2005).

3. UNIQUENESS IN 2-D GRAVITATIONAL INVERSION

Now specialize the above 3-D formulation to a 2-D shell of infinitesimal thickness of arbitrary shape with surface density distribution σ(r). The majority of the mass transports of interest to GRACE TVG studies do occur on the surface of the Earth. The multipole (4) now reduces to the surface integral over the unit sphere:

along with its Eulerian temporal variation:

where R(Ω) is the radius along any radial direction Ω as prescribed by the shape (‘topography’) of the thin shell. The latter is assumed to be singly connected, globally closed surface so that R(Ω) is well-defined, but can always be tailored to stand for regional contributions simply by setting off-region σ to zero. Relation (5) stays valid whether 3-D or 2-D. Strictly speaking, here σ(Ω) designates the density per unit nominal latitude–longitude ‘grid’, rather than per unit area which is actually infeasible to define on a general non-spherical body.

Eq. (8), simple as is, is the key equation for what follows. In practical Earth applications, a complication befalls to it when considering the non-rigidity of the solid Earth which would yield elastically in proportion under the surface load Δσ. Yet the ‘apparent’ Δσ(Ω) (e.g. seasonal snow accumulation or melting away of glacial ice by direct observation) typically takes no account of the ensuing loading change in R(Ω) nor the resultant mass redistribution inside the solid Earth that is actually non-surficial. A compound effect, that contribution to the multipole can be represented by proportionality parameters depending on the 3-D structure of the planetary body. In the case of a simplified spherically symmetric Earth model (say PREM), this loading effect is decomposed conveniently into a series of negative fractional load Love numbers |${k^{\prime}}_n$| for each SH degree n (Farrell 1972). In this paper, we shall adopt these idealized load Love numbers for the real Earth to a good approximation. The net outcome onto the multipoles (from ‘apparent’ to ‘true’) is then a multiplicative factor |$(1 + {k^{\prime}}_n)$|⁠, about 0.70 for n = 2 for the Earth and monotonically approaching 1 as n increases. This practice of course has prevailed in all relevant formulae in the literature.

The multipole ΔQ_nm in eq. (8) is again in the form of inner product, this time the 2-D 〈Δσ, Rⁿ⁺²|${\tilde{Y}}_{nm}$|〉. Now the basis functions {Rⁿ⁺²|${\tilde{Y}}_{nm}$|} for ΔQ_nm do form a complete set that spans the Hilbert space in 2-D, because (i) the SH {Y_nm} do so themselves already; (ii) the empirical radius function R(Ω) of an arbitrary body is in general linearly independent of the mathematical functions Y_nm(Ω) and as such would not diminish the latter's degree of freedom. Therefore, in light of Appendix A, eq. (8) (and of course eq. 7 as well) asserts a profound fact: the gravitational inversion in 2-D is unique. In effect, the multipoles collectively contain sufficient information for reconstructing the surface density Δσ(Ω) uniquely, a generalization of a conclusion of Chao (2005) for the special case of spherical surface (see below). For a more general discussion of uniqueness of 2-D inversion in potential theory, see, for exmple, Král (1980). We add that here the SH are utilized as a mathematical tool, irrelevant to the actual practice of GRACE Level-2 data provided in the SH format.

Unlike {Y_nm} itself, the basis {Rⁿ⁺²|${\tilde{Y}}_{nm}$|}, although complete, is not an orthogonal set. That necessitates the inverting of the full-ranked transfer matrix M (see Appendix A) to solve (uniquely) the surface density Δσ(Ω), where the role of {e_k} in Appendix A is played by {Y_nm}. The Δσ(Ω) can then be placed at due locations at the true radius R(Ω) where it belongs.

In GRACE applications, one specifies a certain maximum degree N for the SH expansions in view of the spatial resolution of interest, with corresponding ‘grid size’ (for instance N = 60 with spatial resolution of ∼300 km), giving M the finite dimension (N + 1)² × (N + 1)². The truncation is justified by the Fact that the truncated subspace, that spanned by SH of degree > N, is orthogonal, hence linearly independent, to the retained subspace spanned by SH of degree ≤ N.

4. SPECIAL CASE: SPHERICAL SURFACE

The argument thus far follows direct logic from eq. (8) that involves no reference to any ‘standard’ spherical model. Conversely, it is nevertheless instructional to re-inspect the special case of a spherical thin shell in the above framework. The first studies of TVG long prior to GRACE launch, not surprisingly, targeted the physical processes occurring on Earth's surface which was set to be spherical for simplicity. To obtain the proper equation one simply lets R(Ω) = a in eq. (8) be a constant which then escapes the integral:

Eq. (9), upon switching to the 4π-normalized SH and converting to Stokes coefficients by relation (5), can be re-written (Wahr et al. 1998; their eq. (6) where ‘a’ means Earth's mean radius instead) as:

which, again, is to be reduced by multiplication of the aforementioned loading factor |$(1 + {k^{\prime}}_n)$| if the Δσ here indicates the apparent observations.

Now further expand Δσ(Ω) in SH, which is always possible as {Y_nm} is complete on 2-D, letting (note the complex conjugation):

Substituting it back to eq. (9) immediately gets (cf. Chao 2005):

by virtual of the SH orthonormality relation (3). This neat, seemingly trivial, relation signifies a one-to-one correspondence between the anomalies of gravity vs. surface density, explicitly ensuring the uniqueness of the gravitational inversion on a sphere—as a special case of the general 2-D shell. It is evident via eq. (5) then that on a sphere the gravity Stokes anomaly of each n is weighted down by the factor (2n + 1) in magnitude relative to the surface density anomaly, that is, a density-to-gravity ‘transfer function’ of 1/(2n + 1). In the spirit of Appendix A, the transfer matrix M to be inverted in this case is diagonal, wherein the disappearance of the between-component couplings owes to the orthonormality of the basis functions {Y_nm} for the inner products of eq. (9).

Eq. (11) itself, likewise switched to 4π-normalized SH, now becomes (Wahr et al. 1998; their eq. 14):

on a sphere, to be augmented by division of the loading factor |$(1 + {k^{\prime}}_n)$| in each term to recover the apparent Δσ(Ω).

Eq. (9), first derived by Chao et al. (1987) basing on Chao & Gross (1987), was employed there in evaluating TVG and Earth rotation variations produced by hydrological mass redistributions on a spherical Earth surface. Eq. (13) represents its reciprocal sum over [n, m], and constitutes the well-cited formula of Wahr et al. (1998; their eq. 14) expressed in units of EWT.

5. REALISTIC NON-SPHERICAL SURFACE

In the course of time, higher precisions in modelling and products beyond the spherical-surface approximation became desirable in GRACE TVG research. Indeed eq. (12) would yield somewhat biased values of Δσ(Ω) to be placed on a fictitious sphere while the actual values are supposed to follow the true surface. So the task here is to reconstruct the surface density Δσ on a truer, non-spherical surface given the observed multipoles, which as we have asserted can be achieved uniquely.

Thus, we revert back to eq. (8) that does factor in the surface shape R(Ω) explicitly. For brevity we shall exclude the load Love number factor in the following, meaning in essence that we are treating the true, not the apparent, Δσ(Ω). The factor can of course be readily reinstated whenever desired. Thus, let the generic surface shape R(Ω) be expressed as the superposition of a main sphere of mean radius R₀ plus a small asphericity h << 1 scaled by R₀:

Eq. (8) then gives:

The first equality, or its equivalent form, has served as basis in previous studies, for example, for layered atmospheric effect (Swenson & Wahr 2002; Boy & Chao 2005) as well as topographic effect (Li et al. 2017; Ditmar 2018; Yang et al. 2022) on TVG.

In the second equality, we have invoked the Taylor-linearization by retaining only the first-order O(h) term in the binomial expansion, which is accurate to O(10⁻²) for the Earth even at n as high as 60, especially when h has typically a red spectrum. In the last equality we have inserted the SH expansion (11) for Δσ as before.

The first, dominant term in the last equality of eq. (15), as it should, gives the spherical-surface approximation of eq. (12), to which eq. (15) reduces of course when setting h = 0. It is the basis for Ditmar's (2018) ‘locally spherical’ approximation when the factor R₀ in front, after treated as an assumed constant, is replaced with actual R at the target location during the inversion. As such, it signifies in essence the upward/downward continuation of a Laplace field, akin to the free-air correction conventionally performed on surface gravity measurements. On the same token, the second term, of first-order in h, is akin to the Bouguer correction in gravity.

In case of the Earth, one can conveniently separate R(Ω) into the reference ellipsoidal figure which has an analytical form, plus the ‘apparent’ surface topography that calls for actual elevation models (see below). Strictly speaking their effects on the multipoles do not add linearly because of the nonlinearity of Rⁿ⁺²(Ω) in eq. (8). However, true to O(h) the two contributions do become additive and hence can be evaluated individually. We expect that by and large the ellipsoidal effect will be an order-of-magnitude larger than the topographic effect judging from their relative heights.

The above analytical development is equivalent to that of Ditmar (2018) of inverting for Δσ on a truer Earth surface. In the following we proceed in parallel to Ditmar (2018) but exploit, as alluded to by Chao (2005), the SH multipole formalism to arrive at formulae that are succinct conceptually and numerically for facilitating further utilities.

5.1 On ellipsoidal figure

In the inversion for Δσ, Li et al. (2017) devised an empirical forward calculation in an iterative scheme to reach a final solution for the ‘ellipsoidal correction’ starting from an initial sphere. They found non-negligible numerical differences; for example, eq. (13) would underestimate the true Δσ in polar regions typically by several per cent because the polar regions have shorter radius than a and hence farther distances from the field point that is the satellite orbit.

The polar oblateness due to the centrifugal force of Earth's spin is the predominant term of the geoid anomaly. It can be parametrized with the geometric ellipticity of the reference ellipsoid that best approximates the Earth figure: f = 1 − b/a = 1/298.2572, where a = 6378 km denotes the equatorial radius (i.e. semi-major axis) and b = 6357 km the polar radius (i.e. semi-minor axis). True to first order of f, O(10⁻²), the axi-symmetric ellipsoidal radius function R(Ω) has the form |${\rm{\ }}{R}_0[ {1 - \ \frac{2}{{3\ }}f\ {P}_{20}( {\cos } )} ]$| (e.g. Stacey & Davis 2008), where R₀ = a (1 + f/3)⁻¹ = 6371 km is the mean radius and P₂₀ = (3 cos²θ −1)/2 = (4π/5)^1/2Y₂₀. Hence, |$h( \Omega ) - \frac{2}{{3\ }}f{P}_{20}( {\cos } )$|⁠, and eq. (15) for the ellipsoidal surface becomes:

Each elements in the second term of eq. (16) that comprises the to-be-inverted transfer matrix M (see Appendix A) is in the form of a surface integration of the product of three SHs. A general expression for the latter, with abbreviated notation Y₁ denoting the [n₁, m₁]th SH, etc., reads (e.g. Arfken & Weber 2005; see also Rotenberg et al. 1959):

where the parenthesis denotes the Wigner 3-j symbol, which is subject to a set of selection rules. The specific integrals in eq. (16) involve |$(_{0}^{2}\,\,\,\,_{0}^{n}\,\,\,\,_{0}^{j}) (_{0}^{2}\,\,\,\,_{m}^{n}\,\,\,\,_{k}^{j})$|⁠, which is non-vanishing only when (i) j is between or equal to n ± 2; (ii) j has the same evenness/oddness as n; and (iii) k = m (which is evident via the integration of sinusoids over λ), so matrix M ends up sparse and block-diagonal. The pertinent 3-j symbols can be evaluated term by term via existing tables or online calculators (e.g. Rotenberg et al. 1959; https://www.volya.net/index.php?id=vc). The ensuing (N + 1)² × (N + 1)² transfer matrix M for this purpose:

is a constant matrix, so is its inversion M⁻¹, that need be done once and for all. That completes the task of constructing exactly all Δσ_jk, and so Δσ itself, on the ellipsoidal surface from the given multipoles ΔQ_nm (or the Stokes coefficients) out to maximum degree N.

Alternatively, Ghobadi-Far et al. (2019; see also Li et al. 2017) took on a conceptually distinct approach, via the linear transformation between ellipsoidal and SH expansions according to Jekeli (1988). They derived a straight one-to-one correspondence of the ellipsoidal harmonic components of Δσ to the Stokes coefficients expressed in the ellipsoidal coordinates, reminiscent of the spherical-surface case of Section 4; the computational loads now bears on dealing with the hypergeometric function in the transformation process.

5.2 On topography

For the purpose of surface topography h(Ω) one needs a land terrain elevation model (many to choose from). On the entire ocean area the ocean surface complies well to the geoid besides small geostrophic dynamic height, departing from the reference ellipsoid by no more than ∼100 m anyway.

To proceed further from eq. (15), one can expand h(Ω) in SH like in eq. (11). The numerical scheme for the ensuing transfer matrix M becomes rather tedious, as each element itself now involves a double summation of Wigner 3-j symbol products in the form of eq. (17). Nevertheless, the inversion principle is the same as above. Alternatively, Ditmar (2018) chose to implement the ‘locally spherical’ approximation, whereas Yang et al. (2022) employed an empirically iterative scheme with scaling factors answering to the upward/downward continuation, to reach satisfactory Δσ(Ω) solutions placed on Earth topography.

Finally, as an interesting digression at this point, we examine a homogeneous 3-D body with uniform density ρ₀ and exterior-boundary shape R(Ω). Upon performing the radial integration over r, the 3-D multipole (4) reduces to a 2-D surface integral reminiscent of eq. (7):

Again, further expansion of R(Ω) in SH yields, owing to the power of n + 3, a tedious nonlinear relation between multipoles and the SH coefficients of R(Ω). Such is the case with Chao & Rubincam (1989) and Rubincam et al. (1995) in forward calculating (truncated at reasonable degrees) the gravity of Phobos and Deimos, Mars’ two moons. If the body happens to be near-spherical with small topography h <<1, one can Taylor-linearize Rⁿ⁺³ as above, and eq. (19) becomes rather simple: besides the monopole Q₀₀ = (4π)^−1/2 M (as always), one gets:

that is, a one-to-one correspondence between gravity versus shape anomalies, reminiscent of the spherical-surface eq. (12) for gravity versus density anomalies, likewise with a transfer function of 1/(2n + 1) signifying a ‘smoothing’ effect between Stokes coefficient and density anomaly. As such, the same reasoning as before establishes the uniqueness of mutual conversion between gravity and shape for a near-spherical, uniform-density, 3-D body.

6. DISCUSSION AND CONCLUDING REMARKS

We use the multipole formalism to revisit the Newtonian gravity produced externally by an arbitrary central body. It leads succinctly to exact, general expressions that unify the relevant equations in the literature of GRACE TVG research. We then address the gravitational inversion problem by casting the multipoles in the framework of Hilbert functional space, and raise contrasts between 3-D body (Earth) versus 2-D thin shell (Earth surface), exact versus special-case approximations, non-uniqueness versus uniqueness, and spherical versus realistic Earth surfaces. We assert the uniqueness of the gravitational inversion on an arbitrary 2-D thin shell, thereby the surface mass (re)distribution on an arbitrary Earth surface can be uniquely determined from observed (changes of) multipoles or Stokes coefficients. We stress here that we are not tasking to assess the gravitational effect of the surface topography, but rather to solve for the surface density from gravity given the topography.

An outcome of this uniqueness has been the establishment of the GRACE surface mascon/EWT solutions that is gaining increasing popularity in the course of time. Conversely, it should be remembered that the said 2-D uniqueness is what justifies the mascon/EWT solutions to be physically meaningful in the first place. A standing commonly practice established in eq. (12) or the equivalent eqs (9) and (10) for years, a spherical Earth sees a one-to-one correspondence providing a comfortable conversion from GRACE Level-2 SH solution to surface mass, often given in units of EWT.

The GRACE surface mascon fields have been released by, primarily, NASA/Goddard Space Flight Center (GSFC), CSR of the University of Texas, and JPL (Tellus). As Level-3 Products, the mascon solutions are derived directly from Level-1B inter-satellite orbital ranging data, bypassing Level-2 TVG solutions in the SH form. The JPL mascon solutions are inverted from orbit data via the partials linking them (Watkins et al. 2015), whereas the GSFC (Rowlands et al. 2010; Luthcke et al. 2013; Loomis et al. 2019) and CSR (Save et al. 2016) mascon solutions result through the chain rule via orbit-to-multipole with multipole-to-mascon partials. Intrinsically all these mascon solutions are performed with respect to a spherical surface, stemming from or equivalent to eqs (9) and (10); so the ensuing products as provided are understood to be a spherical-surface approximations.

Yet to better account for a truer Earth surface, the GSFC solution has further considered the geometrical upward/downward continuation (SH component-wise) by replacing the constant a in eq. (9) with the variable ellipsoidal R of the target location during the inversion (B. Loomis, personal communication, 2022). The CSR solutions have likewise been upgraded recently from sphere to ellipsoid by adopting in post-processing Ditmar's (2018) locally spherical approximation (see http://www2.csr.utexas.edu/grace/RL06_mascons.html). While the above two ellipsoidal fixes are equivalent, we maintain that the scheme proposed in eq. (15) here offers an exact and convenient form for the multipole-to-mascon partials readily implementable in the inversion.

Regardless of the actual form, the establishment of the surface mascon/EWT solutions from GRACE observations naturally raises the caveat (Chao 2016): when using such products one is making the implicit presumption that the mass transports do take place on and only on Earth's surface, spherical or otherwise. In principle the presence of any internal mass redistribution invalidates eq. (8)'s application to the real Earth; indiscriminate usage of mascon can lead to misinterpretation of the geophysical truth.

Finally, coming to bear is the question as to what constitutes a ‘surface’ in reality, if not a mathematical infinitesimal thin shell. Given that the surface density can be destined duly to the true Earth surface, the surface mass-transport processes can comfortably include those taking place in the thin oceans, cryosphere and land hydrology, in fact the majority of GRACE TVG research interests thus far. The atmospheric with its layered structure across thickness of tens of kilometres deserves more attention (Swenson & Wahr 2002; Boy & Chao 2005), something of interest to the de-aliasing pre-processing of GRACE data. In any event the assessment of the precision and error propagation when considering height dependence lies in the factor Rⁿ⁺² according to eq. (8).

Non-surficial large-scale processes abound inside the Earth (cf. Chao 2016). The most prominent are probably the bodily tides (including the pole tide) and the glacial isostatic adjustment; as such, the GRACE mascon/EWT releases are removed of their effects to the extent of the state-of-the-art models. Besides the aforementioned surface loading deformations, earthquakes are also significant non-surficial process as seismic deformations do run deep. Other internal mass-transport phenomena on large scales include plate tectonics-related movements, fluid core convections and inner-core motions; they are characterized typically by slow or secular temporal trends.

Acknowledgement

I thank Jin Li, ShengAn Shih, Bryant Loomis and Pavel Ditmar for valuable and informative discussions in the course of this work, and Khosro Ghobadi-Far for an insightful review. The work is supported by Taiwan Ministry of Science and Technology, under grant no. 111–2116-M-001–024.

DATA AVAILABILITY

This paper is theoretical in nature, and involves no real observational data, and no data sets were generated or analysed.

References

APPENDIX A: ON (NON-)UNIQUENESS OF LINEAR INVERSION

The generic linear inverse problem undern consideration in this paper can be readily cast in the framework of Hilbert function space with the notion of inner product, or ‘projection’. Simply put, the question is: is it possible to construct the whole from a number of projections?

Suppose we are given all the inner products of the (generally complex-valued) function F(x) onto a certain basis function set {B_j(x); j = 1,2,3,…}, that is, ∫ F(x)B_j*(x) dx integrated over the entire x domain (* denoting complex conjugate), or 〈F, B_j〉 in shorthand. On what conditions is the collective information contained in the set {〈F, B_j〉} sufficient to reconstruct, or ‘invert’ for, F(x) itself uniquely?

Suppose that the basis function set {B_j(x)} forms a complete set for the Hilbert space H, i.e. {B_j(x)} spans H in the sense that any function in H can be expressed as a linear combination of B_j(x). Suppose furthermore {B_j(x)} is orthonormal, that is, 〈B_k, B_l〉 = δ_kl, hence all B_j’s are linearly independent. Then, one can readily impose on any F(x) the expansion:

(For manipulation of linear operator and inner product in Hilbert space, see, e.g. Parker 1994; also Chao 1982.) Conversely, it follows that knowing 〈F, B_j〉 for all j gives back F(x) fully by eq. (A1). We say that this inversion is unique. Case in point in the context of this paper is the multipole Q_nm of eq. (9) on a spherical surface, expressed in the form of inner product 〈σ, Y_nm〉 of surface density σ(Ω) (playing the role of F(x)) with Y_nm(Ω) (playing the role of B_j(x)).

The argument above for uniqueness is valid even when the set {B_j(x)}, being complete and linearly independent, is not orthogonal. This can be shown explicitly as follows. One can always construct through the Gram–Schmidt process an orthonormal basis function set, call it {e_k(x)}, by certain linear combinations of the complete {B_j(x)} that spans H. Now expand both F(x) and B_j(x) into components of {e_j(x)} which, like eq. (A1), is always possible given the completeness of the orthonormal {e_k(x)}:

The given projection 〈F, B_j〉, subject to the orthonormality of {e_k(x)}, can now be written as:

The inner products 〈B_j, e_k〉 in eq. (A3) constitute a ‘transfer’ matrix M linearly relating the projections 〈F, B_j〉 with 〈F, e_k〉. M is full-ranked because both {B_j(x)} and {e_k(x)} are linearly independent. When inverted into M⁻¹, that yields the unique solution of the components 〈F, e_k〉, and hence F(x) itself, via eq. (A2). Case in point is the multipole Q_nm of eq. (7) in the form of inner product 〈σ, Rⁿ⁺²Y_nm〉 on a 2-D thin shell.

If, on the other hand, {B_j(x)} does not form a complete set that spans H, then the whole argument above breaks down. For one thing, eq. (A1) would not even exist. Even though eq. (A2) is always allowed, the transfer matrix M formed of such 〈B_j, e_k〉 is no longer full-ranked and hence, with zero determinant, cannot be inverted. The problem is under-determined, with insufficient independent information for a full reconstruction of F(x). Another way to state this is in terms of non-uniqueness for the inversion: The set of projections 〈F, B_j〉 collectively only contains enough information to reconstruct a partial F(x)—the projection of F(x) onto the proper sub-space of H that {B_j(x)} do span. As such, any function that happens to have the same subspace projection but different in the null subspace, wherein all elements are orthogonal to {B_j(x)}, is a legitimate solution, but not a unique solution. We say that this inversion is non-unique. Case in point is the multipole Q_nm of eq. (4) in the form of inner product 〈ρ, rⁿY_nm〉 in 3-D, for which the solid harmonics {rⁿY_nm} are known to be incomplete (see Section 2.3).

The Hilbert space H being infinite-dimensional, all summations above in principle range up to index order infinity, which of course is unrealistic. In practice one is allowed a truncation of the summations out to certain finite index order K, limited typically by the spatial or temporal resolution as the case may be. That makes matrix M finite-dimensional at K × K within the K-dimensional proper subspace of H, and invertible numerically, unless rank-deficient.

While we have used single-variable scalar functions for illustration above, the same reasoning holds for multivariable, vector/tensor functions where the inner product is defined accordingly, and regardless of the physical dimension of the Euclidean space in which it works, whether 3-D or 2-D or otherwise.