The Accuracy Assessment of Lithospheric Density Models

1. Introduction

Gravimetric forward and inverse modelling techniques are essential numerical tools applied in physical geodesy and gravimetric geophysics. In physical geodesy applications, these methods are used to compute topographic and terrain gravity corrections in a gravimetric geoid modelling—e.g., [1,2,3,4,5,6,7,8,9,10,11] and compile isostatic gravity maps—e.g., [12]. In gravimetric geophysics, these methods are used to compile Bouguer and mantle gravity maps—e.g., [13,14,15,16,17,18]. Furthermore, numerous techniques have been developed and applied for a gravimetric interpretation of the Earth’s inner structure—e.g., [19,20,21].

Whereas the gravimetric inverse modelling is applied to determine an unknown density structure or density interface from observed gravity data, the gravimetric forward modelling is used to compute gravitational field quantities generated by a known density structure or density interface. Several different methods have been developed for a local and regional gravimetric forward modelling based on solving the Newton’s volume integral in the spatial domain by applying numerical, semi-analytical, and analytical methods—e.g., [2,[22,23,24,25,26,27,28,29,30,31,32]]. In global (and large-scale regional) applications, methods of solving the Newton’s volume integral in the spectral domain utilize a spherical harmonic analysis and synthesis of gravity and density structure models—e.g., [33,34,35,36,37,38,39,40,41,42].

The validation of accuracy and numerical efficiency of newly developed methods for a gravimetric forward and inverse modelling is often done by comparing results with solutions obtained from existing (and typically well-established) numerical methods and procedures. Alternatively, simple or more refined synthetic density models (designed for various configurations of different geometrical density bodies) have been more recently used for testing and validation of spatial methods developed for a local and regional gravimetric forward modelling and inversion by adopting a planar approximation. For global (or large-scale regional) applications, synthetic density models should mimic more realistically the Earth’s actual shape and inner density structure. Such synthetic density models have already been used to validate numerical techniques involved in a gravimetric geoid modelling—e.g., [43,44,45]. Other examples of possible applications could be given in studies of the sediment bedrock morphology—e.g., [46], the lithospheric and mantle density structure—e.g., [16,47], the crustal thickness—e.g., [16,18,[47,48,49,50,51,52], the dynamic and residual topography—e.g., [53,54,55], or the oceanic lithosphere thermal contraction and its isostatic rebalance—e.g., [56]. To construct a global synthetic density model that closely resembles the Earth’s shape and inner structure, available global topographic and density structure models could be used for this purpose together with additional models that provide information about the Earth’s inner structure (such as crust and lithospheric thickness models). Moreover, constraints have to be applied so that the mass of Earth’s synthetic model is equal to the Earth’s total mass (excluding the atmosphere), and the gravitational field generated by the Earth’s synthetic model closely agrees with the Earth’s gravitational field (in terms of a geoidal geometry, gravity, and gravity gradient).

A number of seismic velocities and mass density models have been developed based on the analysis of tomographic data, while incorporating geophysical, geochemical, and geothermal constraints—e.g., [57,58,59,60,61,62,63,64]; for an overview of these models see also Trabant [65]. A practical application of global 1-D density models, such as the PREM [58] or AK135-F [60], is limited by the absence of a lateral density information. To address this issue, 1-D reference density models could be refined by incorporating 2-D or 3-D global lithospheric and mantle density models to achieve a more realistic representation of the Earth’s inner density structure. Whereas reliable 3-D mantle density models are rare, a number of 3-D crustal and lithospheric density models have been developed and published. Nataf and Ricard [66] derived the crustal and upper-mantle density model based on the analysis of seismic data and additional constrains such as heat flow and chemical composition. Mooney et al. [67] compiled the CRUST5.0 global crustal model with a 5°×5° spatial resolution. Later, the updated global crustal model CRUST2.0 was compiled with a 2°×2° resolution by Bassin et al. [68]. The CRUST1.0 is the most recent version, complied globally with a 1°×1° resolution—e.g., [69]. The CRUST2.0 and CRUST1.0 incorporate also a lateral density structure within the uppermost mantle. Pasyanos et al. [70] compiled the LITHO1.0 global lithospheric model (including the asthenosphere). This model was prepared to fit the high-resolution (Love and Rayleigh) surface wave dispersion maps by using the CRUST1.0 crust data and the LLNL-G3D upper mantle model [71] as the a priori information. Compared to similar 3-D density and velocity models, this model provides also information on the lithosphere-asthenosphere boundary (LAB). Hirt and Rexer [72] constructed the Earth2014 global model consisting of topographic, bathymetric, inland bathymetric, and polar glacier bedrock relief datasets. Chen and Tenzer [73] compiled the Earth’s Spectral Crustal Model 180 (ESCM180) by augmenting the Earth2014 and CRUST1.0 models.

Since many parts of the world are not yet sufficiently covered by tomographic surveys, the refinement of 1-D reference density models by incorporating 2-D or 3-D global lithospheric and mantle density models is not simple. Moreover, the direct relation between seismic velocities and mass densities does not exist because a density distribution depends on many other factors (such us temperature, mineral composition, and pressure). In spite of these practical and theoretical limitations, the development of synthetic density models that more or less realistically approximate the Earth’s shape and inner density structure is essential for a more comprehensive assessment of gravimetric forward and inverse modelling techniques (then that based on using simplistic, typically geometric density models). Moreover, input data uncertainties should be known in order to realistically assess the accuracy of synthetic density models.

To inspect possibilities of refining the Earth’s synthetic density model, the accuracy of published lithospheric density models is assessed based on a novel approach presented in this study. This novel approach allows to formulate the error propagation and consequently provide rough error estimates in a gravimetric forward modelling that are attributed to lithospheric density and geometry uncertainties. A gravimetric forward modelling is applied to determine a long-wavelength geoidal geometry that is corrected for gravitational contributions of lithospheric density and thickness variations. Expressions for a gravimetric forward modelling are presented in Section 2, and then used to derive the error analysis in Section 3. Numerical results are briefly presented in Section 4, and the accuracy of results is discussed in Section 5. Major findings are concluded in Section 6.

2. Numerical model

We applied methods for a spherical harmonic analysis and synthesis of gravitational and lithospheric density models to compute the sub-lithospheric mantle geoid. Details on theoretical aspects are given below.

2.1. Geoid

The geoid height $N$ is defined by—e.g., [74]

$$N(\Omega )=\frac{T\left(r_g,\Omega \right)}{{\gamma }_0\left(\phi \right)}$$

(1)

where the disturbing potential $T$ (i.e., the difference between values of the actual and normal gravity potential $W$ and $U$ respectively; $T=W-U$) is stipulated at the geoid surface $\left(r_g,\Omega \right)$. The normal gravity ${\gamma }_0$ in Eq. (1) is computed at the ellipsoid surface according to Somigliana-Pizzetti’s theory [75,76] for the GRS80 [77] parameters. The 3-D position in Eq. (1) and thereafter is defined in the spherical coordinate system $\left(r,\Omega \right)$, where $r$ is the radius and $\Omega =\left(\phi ,\lambda \right)$ denotes the spherical direction with the spherical latitude $\phi$ and longitude $\lambda$. In the context of interpreting long-wavelength features in the geoidal geometry attributed to a sub-lithospheric mantle density structure, the disturbing potential $T$ in Eq. (1) can approximately (cf. Wieczorek [19,78]; Tenzer et al., [79]) be computed by using the following expression—e.g., [74]

$$T\left(r_g,\Omega \right)\cong \frac{\mathrm{GM}}{\mathrm{R}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{m=-n}^n{\mathrm{T}}_{\mathrm{n},\mathrm{m}}{Y}_{\mathrm{n},\mathrm{m}}(\Omega )}}$$

(2)

where $\mathrm{GM}=3986005\times {10}^8$m³ s^-2 is the geocentric gravitational constant, $\mathrm{R}=6371\times {10}^3$ m is the Earth’s mean radius, ${\mathrm{T}}_{\mathrm{n},\mathrm{m}}$ are the (fully-normalized) disturbing potential coefficients, $Y_{\mathrm{n},\mathrm{m}}$ are the (fully-normalized) surface spherical functions of degree n and order m, and $\overline{n}$ is the upper summation index of spherical harmonics.

2.2. Sub-lithospheric mantle geoid

In the forward modelling scheme, the gravitational contribution of lithospheric density heterogeneities is subtracted from a geoidal geometry in order to enhance a gravitational signature of density structures within the mantle below the lithosphere (i.e., within the sub-lithospheric mantle). This procedure yields the sub-lithospheric mantle geoid $N^{\mathrm{SM}}$, defined by

$$N^{\mathrm{SM}}(\Omega )=N(\Omega )+\frac{T^{\mathrm{SM}}\left(r_g,\Omega \right)}{{\gamma }_0\left(\varphi \right)}$$

(3)

The computation of the sub-lithospheric mantle disturbing potential $T^{\mathrm{SM}}$ in Eq. (3) is realized in numerical steps explained next.

2.2.1. Bouguer disturbing potential

The Bouguer disturbing potential $T^{\mathrm{B}}$ is obtained from the disturbing potential $T$ after subtracting the gravitational potentials of topography $V^T$, bathymetry $V^B$, and polar glaciers $V^I$. We then write (cf. Tenzer et al., [48,50])

$$T^{\mathrm{B}}=T-V^T-V^B-V^I$$

(4)

The gravitational potentials of atmosphere [80] and lakes [81] are completely negligible when compared to values of the gravitational potentials $V^T$, $V^B$, and $V^I$.

The topographic potential in Eq. (4) is often computed individually for a uniform and anomalous topographic density. The topographic potential $V^{T,{\rho }^T}$ for a uniform topographic density reads [12]

$$V^{T,{\rho }^{\mathrm{T}}}\left(r_g,\Omega \right)=\frac{\mathrm{GM}}{\mathrm{R}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{m=-n}^n{\mathrm{V}}_{\mathrm{n},\mathrm{m}}^{{\mathrm{T},\mathsf{\rho }}^{\mathrm{T}}}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}}$$

(5)

The potential coefficients ${\mathrm{V}}_{\mathrm{n},\mathrm{m}}^{{\mathrm{T},\mathsf{\rho }}^{\mathrm{T}}}$ in Eq. (5) are given by

$${\mathrm{V}}_{\mathrm{n},\mathrm{m}}^{{\mathrm{T},\mathsf{\rho }}^{\mathrm{T}}}=\frac{3}{2n+1}\frac{{\mathsf{\rho }}^{\mathrm{T}}}{{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{1}{k+1}\frac{{\mathrm{H}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{k}+1\right)}}{{\mathrm{R}}^{k+1}}$$

(6)

where ${\overline{\mathsf{\rho }}}^{\mathrm{Earth}}=5500$ kg.m^-3 is the Earth’s mean mass density, and ${\mathsf{\rho }}^{\mathrm{T}}$ is the average (constant) topographic density. The Laplace harmonics ${\mathrm{H}}_{\mathrm{n}}$ of topographic heights $H$ are defined by the following integral convolution

$$H(\Omega )={\displaystyle \sum _{n=0}^{\infty }{\mathrm{H}}_{\mathrm{n}}(\Omega )},$$

$${\mathrm{H}}_{\mathrm{n}}(\Omega )=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }{\mathsf{\rho }}^{\mathrm{T}}{H}^{\prime }{\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime }={\displaystyle \sum _{m=-n}^n{\mathrm{H}}_{\mathrm{n},\mathrm{m}}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )},$$

(7)

where ${\mathrm{H}}_{\mathrm{n},\mathrm{m}}$ are the topographic coefficients, ${\mathrm{P}}_{\mathrm{n}}$ is the Legendre polynomial for the argument $t$ of cosine of the spherical angle $\psi$ between two points $\left(r,\Omega \right)$and $\left(r^{\prime },{\Omega }^{\prime }\right)$; i.e., $t=\cos \psi$. The infinitesimal surface element on the unit sphere is denoted as $d{\Omega }^{\prime }=\cos {\phi }^{\prime }d{\phi }^{\prime }d{\lambda }^{\prime }$, and $\Phi =\left\{{\Omega }^{\prime }=\left({\phi }^{\prime },{\lambda }^{\prime }\right):{\phi }^{\prime }\in \left[-\mathsf{\pi }/2,\mathsf{\pi }/2\right]\wedge {\lambda }^{\prime }\in [0,2\mathsf{\pi })\right\}$ is the full spatial angle. The corresponding higher-order terms {${\mathrm{H}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{k}\right)}:k=2,3,\mathrm{...}$} read

$${\mathrm{H}}_{\mathrm{n}}^{\left(\mathrm{k}\right)}(\Omega )=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }{H^{\prime }}^k{\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime }={\displaystyle \sum _{m=-n}^n{\mathrm{H}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{k}\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}$$

(8)

The topographic potential $V_g^{T,\delta {\rho }^T}$ for an anomalous lateral topographic density is computed from

$$V_{}^{T,\delta {\rho }^T}\left(r_g,\Omega \right)=\frac{\mathrm{GM}}{\mathrm{R}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{m=-n}^n{\mathrm{V}}_{\mathrm{n},\mathrm{m}}^{T,\delta {\rho }^T}{Y}_{\mathrm{n},\mathrm{m}}(\Omega )}}$$

(9)

where

$${\mathrm{V}}_{\mathrm{n},\mathrm{m}}^{T,\delta {\rho }^T}=\frac{3}{2n+1}\frac{1}{{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{1}{k+1}\frac{{\mathsf{\delta }\mathsf{\rho }}^{\mathrm{T}}{\mathrm{K}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{k}+1\right)}}{{\mathrm{R}}^{\mathrm{k}+1}}$$

(10)

The topographic density-heights coefficients {${\mathsf{\delta }\mathsf{\rho }}^{\mathrm{T}}{\mathrm{K}}_{\mathrm{n}}^{\left(\mathrm{k}\right)}:k=1,2,3,\mathrm{...}$} in Eq. (10) are given by

$${\mathsf{\delta }\mathsf{\rho }}^{\mathrm{T}}{\mathrm{K}}_{\mathrm{n}}^{\left(\mathrm{k}\right)}(\Omega )=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }\delta {\rho }^T\left({\Omega }^{\prime }\right){K^{\prime }}^k{\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime }={\displaystyle \sum _{m=-n}^n{\mathsf{\delta }\mathsf{\rho }}^{\mathrm{T}}{\mathrm{K}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{k}\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}$$

(11)

The height $K$ in Eq. (11) equals the topographic height $H$ except for areas covered by polar glaciers. In this case, the height $K$ is defined as the topographic height minus the ice thickness. The anomalous lateral topographic density $\delta {\rho }^T$ in Eq. (11) is taken with respect to the average topographic density ${\mathsf{\rho }}^{\mathrm{T}}$, so that $\delta {\rho }^T(\Omega )={\mathsf{\rho }}^{\mathrm{T}}-{\rho }^T(\Omega )$.

The bathymetric potential $V^B$ in Eq. (4) is defined by [82,83]

$$V^B\left(r_g,\Omega \right)=\frac{\mathrm{GM}}{\mathrm{R}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{m=-n}^n{V}_{\mathrm{n},\mathrm{m}}^{\mathrm{b}}{Y}_{\mathrm{n},\mathrm{m}}(\Omega )}}$$

(12)

The bathymetric coefficients $V_{\mathrm{n},\mathrm{m}}^{\mathrm{b}}$ read

$$V_{\mathrm{n},\mathrm{m}}^{\mathrm{b}}=\frac{3}{2n+1}\frac{{\mathsf{\Delta }\mathsf{\rho }}_0^{\mathrm{w}}}{{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}\left(F_{\mathrm{n},\mathrm{m}}^{\mathrm{b}}-\frac{a_1\beta }{{\mathsf{\Delta }\mathsf{\rho }}_0^{\mathrm{w}}}\mathrm{R}{F}_{\mathrm{n},\mathrm{m}}^{\mathrm{b}1}-\frac{a_2\beta }{{\mathsf{\Delta }\mathsf{\rho }}_0^{\mathrm{w}}}{\mathrm{R}}^2{F}_{\mathrm{n},\mathrm{m}}^{\mathrm{b}2}\right)$$

(13)

where ${\mathsf{\Delta }\mathsf{\rho }}_0^{\mathrm{w}}$ is the nominal value of the surface seawater density contrast; and $\beta$, $a_1$, and $a_2$ denote parameters of the depth-dependent seawater density developed by Gladkikh and Tenzer [84]. The numerical coefficients $F_{\mathrm{n},\mathrm{m}}^{\mathrm{b}}$ , $F_{\mathrm{n},\mathrm{m}}^{\mathrm{b}1}$ , and $F_{\mathrm{n},\mathrm{m}}^{\mathrm{b}2}$ in Eq. (13) are computed from

$$F_{\mathrm{n},\mathrm{m}}^{\mathrm{b}}\cong \frac{D_{\mathrm{n},\mathrm{m}}}{\mathrm{R}}-\left(n+2\right)\frac{D_{\mathrm{n},\mathrm{m}}^{\left(2\right)}}{2{\mathrm{R}}^2}+\left(n+2\right)\left(n+1\right)\frac{D_{\mathrm{n},\mathrm{m}}^{\left(3\right)}}{6{\mathrm{R}}^3}$$

$$-\left(n+2\right)\left(n+1\right)n\frac{D_{\mathrm{n},\mathrm{m}}^{\left(4\right)}}{24{\mathrm{R}}^4}+\left(n+2\right)\left(n+1\right)n\left(n-1\right)\frac{D_{\mathrm{n},\mathrm{m}}^{\left(5\right)}}{120{\mathrm{R}}^5}$$

$$+{\displaystyle \sum _{k=5}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{k+1}\frac{D_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}}{{\mathrm{R}}^{\mathrm{k}+1}} ,$$

(14)

$$F_{\mathrm{n},\mathrm{m}}^{\mathrm{b}1}\cong \frac{D_{\mathrm{n},\mathrm{m}}^{\left(2\right)}}{2{\mathrm{R}}^2}-\left(n+2\right)\frac{D_{\mathrm{n},\mathrm{m}}^{\left(3\right)}}{3{\mathrm{R}}^3}+\left(n+2\right)\left(n+1\right)\frac{D_{\mathrm{n},\mathrm{m}}^{\left(4\right)}}{{8\mathrm{R}}^4}$$

$$-\left(n+2\right)\left(n+1\right)n\frac{D_{\mathrm{n},\mathrm{m}}^{\left(5\right)}}{{30\mathrm{R}}^5}+{\displaystyle \sum _{k=4}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{k+2}\frac{D_{\mathrm{n},\mathrm{m}}^{\left(k+2\right)}}{{\mathrm{R}}^{\mathrm{k}+2}}$$

(15)

$$+{\displaystyle \sum _{k=3}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{k+3}\frac{D_{\mathrm{n},\mathrm{m}}^{\left(k+3\right)}}{{\mathrm{R}}^{\mathrm{k}+3}}$$

(16)

The coefficients $D_{\mathrm{n},\mathrm{m}}$ of global bathymetric model are defined by

$${\mathrm{D}}_{\mathrm{n}}(\Omega )=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }D\left({\Omega }^{\prime }\right)P_n\left(\cos \psi \right)d{\Omega }^{\prime }}={\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{\mathrm{D}}_{\mathrm{n},\mathrm{m}}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )},$$

(17)

and their higher-order terms $\left\{{\mathrm{D}}_{\mathrm{n}}^{\left(\mathrm{i}\right)}:i=2,3,4\right\}$ read (ibid.)

$${\mathrm{D}}_{\mathrm{n}}^{\left(\mathrm{i}\right)}(\Omega )=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }{D}^k\left({\Omega }^{\prime }\right)P_n\left(\cos \psi \right)d{\Omega }^{\prime }}={\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{\mathrm{D}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{i}\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )},$$

(18)

where $D$ denotes the bathymetric depth.

The ice potential $V^I$ in Eq. (4) is computed from (cf. Foroughi and Tenzer [85])

$$V_{}^I\left(r_g,\Omega \right)=\frac{\mathrm{GM}}{\mathrm{R}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{m=-n}^n{\mathrm{V}}_{\mathrm{n},\mathrm{m}}^{\mathrm{I}}{Y}_{\mathrm{n},\mathrm{m}}(\Omega )}}$$

(19)

The ice coefficients ${\mathrm{V}}_{\mathrm{n},\mathrm{m}}^{\mathrm{I}}$ in Eq. (19) read

$$V_{\mathrm{n},\mathrm{m}}^{\mathrm{I}}=\frac{3}{2n+1}\frac{{\mathsf{\delta }\mathsf{\rho }}^{\mathrm{I}}}{{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}\left(F{u}_{\mathrm{n},\mathrm{m}}^{}-F{l}_{\mathrm{n},\mathrm{m}}^{}\right),$$

(20)

where the numerical coefficients $F{l}_{\mathrm{n},\mathrm{m}}^{}$ and $F{u}_{\mathrm{n},\mathrm{m}}^{}$are given by

$$F{l}_{\mathrm{n},\mathrm{m}}^{}={\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{1}{k+1}\frac{L_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}}{{\mathrm{R}}^{\mathrm{k}+1}}$$

$$F{u}_{\mathrm{n},\mathrm{m}}^{}={\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{1}{k+1}\frac{{\mathrm{H}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}}{{\mathrm{R}}^{\mathrm{k}+1}}$$

(21)

The coefficients $F{l}_{\mathrm{n},\mathrm{m}}^{}$ utilize the spherical lower-bound functions ${\mathrm{L}}_{\mathrm{n}}$ of a volumetric mass density contrast layer and their higher-order terms [82,83]

$$L_{\mathrm{n}}^{\left(k+1\right)}(\Omega )=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }{H}_L^{k+1}\left({\Omega }^{\prime }\right)P_n\left(t\right)d{\Omega }^{\prime }}={\displaystyle \sum _{m=-n}^n{L}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}$$

(22)

Since the upper-bound of glaciers is identical with the topographic surface, the numerical coefficients $F{u}_{\mathrm{n},\mathrm{m}}^{}$ are generated directly from the height coefficients {${\mathrm{H}}_{\mathrm{n}}^{\left(\mathrm{k}\right)}:k=1,2,\mathrm{...}$}. The ice density contrast ${\mathsf{\delta }\mathsf{\rho }}^{\mathrm{I}}$ in Eq. (20) is taken with respect to a uniform topographic density ${\mathsf{\rho }}^{\mathrm{T}}$, i.e. ${\mathsf{\delta }\mathsf{\rho }}^{\mathrm{I}}={\mathsf{\rho }}^{\mathrm{T}}-{\mathsf{\rho }}^{\mathrm{I}}$, where ${\mathsf{\rho }}^{\mathrm{I}}$ denotes the density of polar glaciers.

2.2.2. Crust-stripped disturbing potential

The gravitational potential of consolidated-crust density contrast $V^C$ and the gravitational potential of sediment density contrast $V^S$ are subtracted from the Bouguer disturbing potential $T^{\mathrm{B}}$ in order to remove the gravitational contribution of density heterogeneities within the whole crust. This procedure yields the crust-stripped disturbing potential $T^{\mathrm{CS}}$ (cf. Tenzer et al. [48])

$$T^{\mathrm{CS}}=T^{\mathrm{B}}-V^S-V^C$$

(23)

The potentials $V^S$ and $V^C$ in Eq. (23) are computed according to a method developed by Tenzer et al. [82] that utilizes the expression for a gravitational potential $V$ generated by an arbitrary volumetric mass layer with a variable depth and thickness while having laterally distributed vertical mass density variations. It reads

$$V\left(r_g,\Omega \right)=\frac{\mathrm{GM}}{\mathrm{R}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{m=-n}^n{\mathrm{V}}_{\mathrm{n},\mathrm{m}}{Y}_{\mathrm{n},\mathrm{m}}(\Omega )}}$$

(24)

The potential coefficients $V_{\mathrm{n},\mathrm{m}}$ in Eq. (24) are defined by

$$V_{\mathrm{n},\mathrm{m}}=\frac{3}{2n+1}\frac{1}{{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{i=0}^I}\left(F{l}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{i}\right)}-F{u}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{i}\right)}\right)$$

(25)

where the coefficients {$F{l}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{i}\right)},F{u}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{i}\right)}:i=0,1,\mathrm{...},I$} are computed from

$$F{l}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{i}\right)}={\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{k+1+i}\frac{L_{\mathrm{n},\mathrm{m}}^{\left(k+1+i\right)}}{{\mathrm{R}}^{\mathrm{k}+1}},$$

$$F{u}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{i}\right)}={\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{k+1+i}\frac{U_{\mathrm{n},\mathrm{m}}^{\left(k+1+i\right)}}{{\mathrm{R}}^{\mathrm{k}+1}}$$

(26)

The terms ${\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{\mathrm{L}}_{\mathrm{n},\mathrm{m}}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}}$ and ${\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{\mathrm{U}}_{\mathrm{n},\mathrm{m}}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}}$ in Eq. (26) define the spherical lower-bound and upper-bound laterally distributed radial density variation functions ${\mathrm{L}}_{\mathrm{n}}$ and ${\mathrm{U}}_{\mathrm{n}}$ of degree $n$. These numerical coefficients combine information on the geometry and mass density (or density contrast) distribution within a volumetric layer. The computation of these coefficients is realized to a certain degree of spherical harmonics from discrete data of the spatial mass density distribution (typically provided by means of density, depth and thickness data) of a particular structural component of the Earth’s interior.

The spherical functions ${\mathrm{L}}_{\mathrm{n}}$ and ${\mathrm{U}}_{\mathrm{n}}$ including their higher-order terms {$L_{\mathrm{n}}^{\left(\mathrm{k}+1+\mathrm{i}\right)},U_{\mathrm{n}}^{\left(\mathrm{k}+1+\mathrm{i}\right)}:k=0,1,\mathrm{...};i=1,2,\mathrm{...},I$} in Eq. (26) are defined by

$$L_{\mathrm{n}}^{\left(k+1+i\right)}(\Omega )=\{\begin{array}{c}\begin{array}{l}\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }\rho \left(D_U,{\Omega }^{\prime }\right)D_L^{k+1}\left({\Omega }^{\prime }\right)P_n\left(t\right)d{\Omega }^{\prime }}\\ ={\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{L}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=0\hspace{1em}\hspace{1em}\end{array}\\ \begin{array}{l}\\ \end{array}\\ \begin{array}{l}\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }\beta \left({\Omega }^{\prime }\right){\alpha }_i\left({\Omega }^{\prime }\right)D_L^{k+1+i}\left({\Omega }^{\prime }\right)P_n\left(t\right)d{\Omega }^{\prime }}\\ ={\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{L}_{\mathrm{n},\mathrm{m}}^{\left(k+1+i\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=1,2,\mathrm{...},I\end{array}\end{array}$$

(27)

and

$$U_{\mathrm{n}}^{\left(k+1+i\right)}(\Omega )=\{\begin{array}{c}\begin{array}{l}\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }\rho \left(D_U,{\Omega }^{\prime }\right)D_U^{k+1}\left({\Omega }^{\prime }\right)P_n\left(t\right)d{\Omega }^{\prime }}\\ ={\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{U}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=0\hspace{1em}\hspace{1em}\end{array}\\ \begin{array}{l}\\ \end{array}\\ \begin{array}{l}\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }\beta \left({\Omega }^{\prime }\right){\alpha }_i\left({\Omega }^{\prime }\right)D_U^{k+1+i}\left({\Omega }^{\prime }\right)P_n\left(t\right)d{\Omega }^{\prime }}\\ ={\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{U}_{\mathrm{n},\mathrm{m}}^{\left(k+1+i\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}i=1,2,\mathrm{...},I\end{array}\end{array}$$

(28)

For a specific volumetric layer, the mass density $\rho$ is either constant $\mathsf{\rho }$, laterally-varying $\rho \left({\Omega }^{\prime }\right)$ or - in the most general case - approximated by the laterally distributed radial density variation model by using the following polynomial function (for each lateral column)

$$\rho \left(r^{\prime },{\Omega }^{\prime }\right)=\rho \left(D_U,{\Omega }^{\prime }\right)+\beta \left({\Omega }^{\prime }\right)\sum _{i=1}^I{\alpha }_i\left({\Omega }^{\prime }\right){\left(\mathrm{R}-r^{\prime }\right)}^i$$

$$for\hspace{1em}\mathrm{R}-D_U\left({\Omega }^{\prime }\right)\ge r^{\prime }>\mathrm{R}-D_L\left({\Omega }^{\prime }\right)$$

(29)

The nominal value of lateral density $\rho \left(D_U,{\Omega }^{\prime }\right)$ is stipulated at the depth $D_U$. This density distribution describes a radial density variation (in terms of coefficients {${\alpha }_i:i=1,2,\mathrm{...},I$} and $\beta$) within a volumetric mass layer at a location ${\Omega }^{\prime }$. Alternatively, when modeling a gravitational field of anomalous density structures (in this case, the density contrast of sediments and consolidated crust), the density contrast $\Delta \rho$ of volumetric mass layer relative to the reference density ${\mathsf{\rho }}^{\mathrm{c}}$ is defined by

$$\Delta \rho \left(r^{\prime },{\Omega }^{\prime }\right)={\mathsf{\rho }}^{\mathrm{c}}-\rho \left(r^{\prime },{\Omega }^{\prime }\right)$$

$$=\Delta \rho \left(D_U,{\Omega }^{\prime }\right)-\beta \left({\Omega }^{\prime }\right){\displaystyle \sum _{i=1}^I{\alpha }_i\left({\Omega }^{\prime }\right){\left(\mathrm{R}-r^{\prime }\right)}^i},$$

$$for\hspace{1em}\mathrm{R}-D_U\left({\Omega }^{\prime }\right)\ge r^{\prime }>\mathrm{R}-D_L\left({\Omega }^{\prime }\right)$$

(30)

where $\Delta \rho \left(D_U,{\Omega }^{\prime }\right)$ is the nominal value of lateral density contrast stipulated at the depth $D_U$. Here the reference density ${\mathsf{\rho }}^{\mathrm{c}}$ of homogenous crust is used.

2.2.3. Mantle disturbing potential

To reveal a gravitational signature attributed to a mantle density structure, the Moho geometry signature has to be subtracted from the crust-stripped disturbing potential. This procedure yields the mantle disturbing potential. According to Tenzer et al. [50], the mantle disturbing potential $T^{\mathrm{M}}$ is obtained from the crust-stripped disturbing potential $T^{\mathrm{CS}}$ after subtracting the gravitational potential of Moho geometry $V^{M,\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}$. We then write

$$T^{\mathrm{M}}=T^{\mathrm{CS}}-V^{M,\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}$$

(31)

The potential $V^{M,\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}$ (for the average Moho density contrast $\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}$) is defined in the following form

$$V^{M,\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}\left(r_g,\Omega \right)=3\frac{\mathrm{GM}}{{\mathrm{R}\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{F}_{\mathrm{n},\mathrm{m}}^{\Delta \mathsf{\rho }\mathrm{M}}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}},$$

(32)

where the numerical coefficients $F_{\mathrm{n},\mathrm{m}}^{\Delta \mathsf{\rho }\mathrm{M}}$ are given by

$$F_{\mathrm{n},\mathrm{m}}^{\Delta \mathsf{\rho }\mathrm{M}}=-\frac{{\mathsf{\Delta }\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}{\left(2n+1\right)\left(n+3\right)}{\displaystyle \sum _{k=1}^{n+3}\left(\begin{array}{c}n+3\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{{\mathrm{R}}^k}{\hat{\mathrm{M}}}_{\mathrm{n},\mathrm{m}}^{\left(k\right)}$$

(33)

The Moho depth spherical functions ${\hat{\mathrm{M}}}_{\mathrm{n}}$ and their higher-order terms {${\hat{\mathrm{M}}}_{\mathrm{n}}^{\left(k\right)}:k=2,3,\mathrm{...}$} read

$${\hat{\mathrm{M}}}_{\mathrm{n}}^{\left(k\right)}=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }{M}^k\left({\Omega }^{\prime }\right){\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime }$$

$$={\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{\hat{\mathrm{M}}}_{\mathrm{n},\mathrm{m}}^{\left(k\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}$$

(34)

As seen in Eq. (34), the coefficients {${\mathrm{M}}_{\mathrm{n},\mathrm{m}}^{\left(k\right)}:k=1,2,\mathrm{...}$} are generated from values of the Moho depth $M$ . The Moho density contrast in Eq. (33) is defined as the average lithospheric mantle density ${\mathsf{\rho }}^{\mathrm{LM}}$ from which the (constant) reference crustal density ${\mathsf{\rho }}^{\mathrm{c}}$ is subtracted. Hence

$$\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}={\mathsf{\rho }}^{\mathrm{LM}}-{\mathsf{\rho }}^{\mathrm{c}}$$

(35)

2.2.4. Lithosphere-stripped disturbing potential

The lithosphere-stripped disturbing potential $T^{\mathrm{CL}}$ is obtained from the mantle disturbing potential $T^{\mathrm{M}}$ by subtracting the lithospheric-mantle gravitational potential $V^{LM}$, so that

$$T^{\mathrm{CL}}=T^{\mathrm{M}}-V^{LM}$$

(16)

For a lateral density distribution function within the lithospheric mantle, the lithospheric-mantle gravitational potential $V^{LM}$ is defined by

$$V^{LM}\left(r_g,\Omega \right)=3\frac{\mathrm{GM}}{{\mathrm{R}\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{n=0}^{\overline{n}}\frac{1}{2n+1}}$$

$$\times {\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{{\mathrm{R}}^{\mathrm{k}+1}\left(k+1\right)}{\displaystyle \sum _{m=-n}^n\left(L_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}-{\mathrm{M}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}\right)Y_{\mathrm{n},\mathrm{m}}(\Omega )}$$

(37)

where the Moho coefficients ${\mathrm{M}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}$ describe the geometry and lateral density (contrast) distribution at the Moho interface, and the LAB coefficients ${\mathrm{L}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}$ describe the geometry and lateral density (contrast) distribution at the LAB.

The Moho coefficients ${\mathrm{M}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}$ in Eq. (37) are defined in the following form

$${\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{\mathrm{M}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{k}\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }\delta {\rho }^{LM}\left({\Omega }^{\prime }\right)M^k\left({\Omega }^{\prime }\right){\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime }$$

(38)

where $M$ denotes the Moho depth (see Eq. 34). The (lateral) lithospheric mantle density contrast $\delta {\rho }^{LM}$ in Eq. (38) is computed as the difference between the adopted (lithospheric mantle) reference density ${\mathsf{\rho }}^{\mathrm{LM}}$and the lithospheric mantle density ${\rho }^{LM}$; i.e., $\delta {\rho }^{LM}\left({\Omega }^{\prime }\right)={\mathsf{\rho }}^{\mathrm{LM}}-{\rho }^{LM}\left({\Omega }^{\prime }\right)$.

The LAB coefficients ${\mathrm{L}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}$ in Eq. (37) read

$${\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{\mathrm{L}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{k}\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}=\frac{2n+1}{4\pi }{\displaystyle \underset{\Phi }{\iint }\delta {\rho }^{LM}\left({\Omega }^{\prime }\right)L^k\left({\Omega }^{\prime }\right){\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime }$$

(17)

where $L$ is the LAB depth. If we consider a current resolution of global lithospheric density models (1×1 arc-deg or lower) that corresponds (by means of a half wavelength) to a spectral resolution up to the degree 180 of spherical harmonics, the lithospheric-mantle gravitational potential can sufficiently be computed by using binomial series in Eq. (37) up to only the third-order term.

2.2.5. Sub-lithospheric mantle disturbing potential

The computation of the lithosphere-stripped disturbing potential $T^{\mathrm{CL}}$ in Eq. (36) enhanced a signature of the LAB geometry that has to be removed in order to enhance a gravitational signature of the sub-lithospheric mantle density structure. This procedure is realized by stripping the lithosphere with respect to the density contrast between the reference (lithospheric mantle) density and the asthenosphere density in order to obtain the sub-lithospheric mantle disturbing potential. The sub-lithospheric mantle disturbing potential $T^{\mathrm{SM}}$ is then computed as

$$T^{\mathrm{SM}}=T^{\mathrm{CL}}-V^{LAB},$$

(40)

where the potential $V^{LAB}$ of LAB geometry is defined by

$$V^{LAB}\left(r_g,\Omega \right)=3\frac{\mathrm{GM}}{{\mathrm{R}\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{n=0}^{\overline{n}}\frac{1}{2n+1}}$$

$$\times {\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{{\mathrm{R}}^{\mathrm{k}+1}\left(k+1\right)}{\displaystyle \sum _{m=-n}^n{\widehat{L}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}{Y}_{\mathrm{n},\mathrm{m}}(\Omega )}$$

(41)

The coefficients ${\widehat{L}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}$ in Eq. (41) read

$${\displaystyle \sum _{\mathrm{m}=-\mathrm{n}}^{\mathrm{n}}{\hat{\mathrm{L}}}_{\mathrm{n},\mathrm{m}}^{\left(\mathrm{k}\right)}{\mathrm{Y}}_{\mathrm{n},\mathrm{m}}(\Omega )}=\frac{2n+1}{4\pi }\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}{\displaystyle \underset{\Phi }{\iint }{L}^k\left({\Omega }^{\prime }\right){\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime },$$

(42)

where $L$ is the LAB depth (see also Eq. 39), and the LAB density contrast $\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}={\mathsf{\rho }}^{\mathrm{A}}-{\mathsf{\rho }}^{\mathrm{LM}}$ is defined as the difference between the lithospheric mantle density ${\mathsf{\rho }}^{\mathrm{L}}$ and the asthenospheric density ${\mathsf{\rho }}^{\mathrm{A}}$. We note that the computation of $V^{LAB}$ can again be restricted up to the third-order terms of a binomial series (as in the case of computing $V^{LM}$ in Eq. 37).

3. Model uncertainties

The numerical model presented in Section 2 was used in this section to derive the error propagation in modelling of the sub-lithospheric mantle geoid. According to results presented in Section 4 (see Table 1), the gravitational contributions of lithospheric mantle and LAB are one order of magnitude larger than the corresponding contributions of topographic and crustal density structures. We, therefore, focused only on the estimation of errors due to lithospheric mantle density and lithospheric thickness uncertainties, while disregarding much smaller (in absolute sense) modelling errors due to crustal density uncertainties.

3.1. Errors due to lithospheric thickness uncertainties

In the sub-lithospheric mantle geoid modelling, the lithospheric thickness information is used to compute the lithospheric-mantle gravitational potential $V^{LM}$ (in Eq. 39) and the LAB potential $V^{LAB}$(in Eq. 41). We then established a relation between the error of sub-lithospheric mantle geoid ${\epsilon }_{N^{\mathrm{SM}}}$ and the error of lithospheric thickness ${\epsilon }_L$ in the following form

$${\epsilon }_{N^{\mathrm{SM}}}=-\frac{1}{{\gamma }_0}\left(\frac{\partial V^{LM}}{\partial L}+\frac{\partial V^{LAB}}{\partial L}\right){\epsilon }_L$$

$$=-\frac{3}{{\gamma }_0}\frac{\mathrm{GM}}{{\mathrm{R}\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{n=0}^{\overline{n}}\frac{1}{2n+1}}{\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{{\mathrm{R}}^{\mathrm{k}+1}\left(k+1\right)}{\displaystyle \sum _{m=-n}^n\left(L_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}+{\widehat{L}}_{\mathrm{n},\mathrm{m}}^{\left(k+1\right)}\right)Y_{\mathrm{n},\mathrm{m}}(\Omega )}$$

$$=-\frac{3}{{\gamma }_{0}4\pi }\frac{\mathrm{GM}}{{\mathrm{R}\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{{\mathrm{R}}^{\mathrm{k}+1}\left(k+1\right)}}$$

$$\times {\displaystyle \underset{\Phi }{\iint }\left[\delta {\rho }^{LM}\left({\Omega }^{\prime }\right)+\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}\right]\left(k+1\right)L^k\left({\Omega }^{\prime }\right){\epsilon }_L\left({\Omega }^{\prime }\right){\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime }$$

(43)

After disregarding terms $k\ge 1$ , the expression in Eq. (43) simplifies to

$${\epsilon }_{N^{\mathrm{SM}}}\cong -\frac{3}{{\gamma }_{0}4\pi }\frac{\mathrm{GM}}{{\mathrm{R}}^2{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \underset{\Phi }{\iint }\left[\delta {\rho }^{LM}\left({\Omega }^{\prime }\right)+{\mathsf{\Delta }\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}\right]{\epsilon }_L\left({\Omega }^{\prime }\right){\displaystyle \sum _{n=0}^{\overline{n}}{\mathrm{P}}_{\mathrm{n}}\left(t\right)}d}{\Omega }^{\prime }$$

(44)

For n = 0, we have

$${\epsilon }_{N^{\mathrm{SM}}}\cong -\frac{3}{{\gamma }_{0}4\pi }\frac{\mathrm{GM}}{{\mathrm{R}}^2{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \underset{\Phi }{\iint }\left[\delta {\rho }^{LM}\left({\Omega }^{\prime }\right)+\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}\right]{\epsilon }_L\left({\Omega }^{\prime }\right)d}{\Omega }^{\prime }$$

(45)

When assuming that the largest contribution to the error ${\epsilon }_{N^{\mathrm{SM}}}$comes from lithospheric thickness uncertainties in the vicinity of computation point (for which the error is estimated), the global integration domain in Eq. (45) could be disregarded. Instead, we reduce a functional relation only for a position of computation point. We then write

$${\epsilon }_{N^{\mathrm{SM}}}\approx -\frac{3}{{\gamma }_{0}4\pi }\frac{\mathrm{GM}}{{\mathrm{R}}^2{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}\left(\delta {\rho }^{LM}+\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}\right){\epsilon }_L$$

(46)

Inserting for ${\gamma }_0\approx {\mathrm{GM}/\mathrm{R}}^2$ in Eq. (46), we arrive at

$${\epsilon }_{N^{\mathrm{SM}}}\approx -\frac{3}{4\pi }\left(\frac{\delta {\rho }^{LM}+\Delta {\rho }^{\mathrm{L}/\mathrm{A}}}{{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}\right){\epsilon }_L$$

(47)

3.2. Errors due to lithospheric mantle density uncertainties

The lithospheric mantle density information is used to compute the Moho density contrast $\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}$ (Eq. 35), the lithospheric mantle density variations $\delta {\rho }^{LM}$ (Eq. 38), and the LAB density contrast $\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}$ (Eq. 42). This involves the computation of the Moho potential $V^{M,\Delta {\mathsf{\rho }}^{c/m}}$ (in Eq. 32), the lithospheric-mantle gravitational potential $V^{LM}$ (in Eq. 39), and the LAB potential $V^{LAB}$(in Eq. 41). The error propagation for the lithospheric mantle density uncertainties is then established in the following form

$${\epsilon }_{N^{\mathrm{SM}}}=-\frac{1}{{\gamma }_0}\left(\frac{\partial V^{M,\Delta {\rho }^{c/m}}}{\partial \Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}+\frac{\partial V^{LM}}{\partial \delta {\rho }^{LM}}{\epsilon }_{\delta {\rho }^{LM}}+\frac{\partial V^{LAB}}{\partial \Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}\right)$$

$$=-\frac{3}{{\gamma }_{0}4\pi }\frac{\mathrm{GM}}{{\mathrm{R}\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \sum _{n=0}^{\overline{n}}{\displaystyle \sum _{k=0}^{n+2}\left(\begin{array}{c}n+2\\ k\end{array}\right)}\frac{{\left(-1\right)}^k}{{\mathrm{R}}^{\mathrm{k}+1}\left(k+1\right)}}$$

$$\times {\displaystyle \underset{\Phi }{\iint }\left\{{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}{M}^{k+1}\left({\Omega }^{\prime }\right)+\left[{\epsilon }_{\delta {\rho }^{LM}}\left({\Omega }^{\prime }\right)L^{k+1}\left({\Omega }^{\prime }\right)-{\epsilon }_{\delta {\rho }^{LM}}\left({\Omega }^{\prime }\right)M^{k+1}\left({\Omega }^{\prime }\right)\right]{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}{L}^{k+1}\left({\Omega }^{\prime }\right)\right\}{\mathrm{P}}_{\mathrm{n}}\left(t\right)d}{\Omega }^{\prime },$$

(48)

where ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}$, ${\epsilon }_{\delta {\rho }^{LM}}$ and ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}$ denote, respectively, the errors in the Moho density contrast $\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}$, the anomalous lithospheric mantle density $\delta {\rho }^{LM}$, and the LAB density contrast $\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}$.

After disregarding terms $k\ge 1$ in Eq. (48), we arrive at

$${\epsilon }_{N^{\mathrm{SM}}}\cong -\frac{3}{{\gamma }_{0}4\pi }\frac{\mathrm{GM}}{{\mathrm{R}}^2{\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \underset{\Phi }{\iint }\{{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}{M}^{k+1}\left({\Omega }^{\prime }\right)+{\epsilon }_{\delta {\rho }^{LM}}\left({\Omega }^{\prime }\right)L^{k+1}\left({\Omega }^{\prime }\right)}$$

$$-{\epsilon }_{\delta {\rho }^{LM}}\left({\Omega }^{\prime }\right)M^{k+1}\left({\Omega }^{\prime }\right)+{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}{L}^{k+1}\left({\Omega }^{\prime }\right)\}{\displaystyle \underset{\Phi }{\iint }{\displaystyle \sum _{n=0}^{\overline{n}}{\mathrm{P}}_{\mathrm{n}}\left(t\right)}d}{\Omega }^{\prime }$$

(49)

Considering that n = 0 and ${\gamma }_0\approx {\mathrm{GM}/\mathrm{R}}^2$, we get

$${\epsilon }_{N^{\mathrm{SM}}}\cong -\frac{3}{4\pi {\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}{\displaystyle \underset{\Phi }{\iint }\left\{{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}{M}^{k+1}\left({\Omega }^{\prime }\right)+{\epsilon }_{\delta {\rho }^{LM}}\left({\Omega }^{\prime }\right)L\left({\Omega }^{\prime }\right)-{\epsilon }_{\delta {\rho }^{LM}}\left({\Omega }^{\prime }\right)M\left({\Omega }^{\prime }\right)+{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}L\left({\Omega }^{\prime }\right)\right\}d}{\Omega }^{\prime }$$

(50)

By analogy with Eq. (46), the expression in Eq. (50) is further simplified to

$${\epsilon }_{N^{\mathrm{SM}}}\approx -\frac{3}{4\pi {\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}\left[{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}M+{\epsilon }_{\delta {\rho }^{LM}}\left(L-M\right)+{\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}L\right]$$

(51)

As seen in Eq. (51), the sub-lithospheric mantle geoid error ${\epsilon }_{N^{\mathrm{SM}}}$ in this case depends on the accuracy of anomalous lithospheric mantle density $\delta {\rho }^{LM}$. In addition, the errors ${\epsilon }_{\Delta {\rho }^{\mathrm{L}/\mathrm{A}}}$ and ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}$ in adopted values of the LAB density contrast $\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}$ and the Moho density contrast $\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}$ affect the accuracy.

4. Results

Tenzer and Chen [81] applied numerical procedures (described in Section 2) to compute the global sub-lithospheric mantle geoid globally with a spectral resolution up to degree 180 of spherical harmonics. They used the EIGEN-6C4 [86] global gravitational model, the Earth2014 topographic, bathymetric, and glacial bedrock relief datasets [72], the UNB_TopoDens global lateral topographic density model [87], the total sediment thickness data for the world's oceans and marginal seas [88], the CRUST1.0 [69] global seismic crustal model updated for the sediment and crustal layers of the Antarctic lithosphere by Baranov et al. [89], and the LITHO1.0 [70] global seismic lithospheric model. They adopted the reference density values 2670 kg.m^-3 (cf. Hinze [90]; Artemjev et al. [91)) for the crust above the geoid, 2900 kg.m^-3 for the crust below the geoid, and 3300 kg.m^-3 for the lithospheric mantle. To further improve the accuracy of forward modelling, they defined the ocean density contrast for a depth-dependent seawater density function [47,84], and applied a density model of marine sediments [92] under marginal seas and oceans.

The gravitational contributions of individual lithospheric density structures are plotted in Figure 1, with the statistical summary of results in Table 1. Modifications of the geoidal geometry after subtracting gravitational contributions of individual lithospheric density structures are presented in Figure 2, with the statistical summary in Table 2.

A long-wavelength geoidal geometry (Figure 2a) reflects mainly a mantle density structure—e.g., [93,94], while a signature of crustal density and geometry variations (including a topographic surplus of large orogenic formations) is less prominent. After removing the gravitational contributions topography and lithospheric density and geometry (i.e., LAB) variations, the resulting (sub-lithospheric mantle) geoidal geometry (see Figure 2f) should (optimally) enhance the gravitational signature of sub-lithospheric mantle.

It is well known fact that a spatial pattern in dynamic topography models—e.g., [54,55],[95,96,97,98] mainly reflects a mantle convection flow, with maxima marking locations of the African and South Pacific superplumes that represent large-scale regions characterized by an elevated topography and shallow ocean-floor depts caused by a low-density upwelling mantle material from the core-mantle boundary. The corresponding minima are explained by a high-density mantle downwelling flow. Tomographic studies—e.g., [99] indicate that superplumes coincide with locations of large low-shear-velocity provinces (LLSVPs) beneath Africa (called Tuzo) and the South Pacific (called Jason). Since the gravitational contributions of lithospheric density and thickness variations have been modelled and consequently removed in our forward modelling procedure, the sub-lithospheric mantle geoid should mainly reflect deep mantle density anomalies. Consequently, the spatial pattern of sub-lithospheric mantle geoid (especially its long-wavelength spectrum up to degree 5 of spherical harmonics; see Figure 3) does not closely agree with a spatial pattern of dynamic topography models (cf. Tenzer and Chen [81]). A comparative study of spatial patterns in the sub-lithospheric mantle geoid and dynamic topography models is thus not meaningful for the assessment of errors in our gravimetric forward modelling.

As seen in Figure 3, the long-wavelength pattern in the sub-lithospheric mantle geoid is characterized by a large positive anomaly in the Central Pacific and a less pronounced anomaly in the Central Atlantic with a prolonged shape across the Atlantic Ocean and further extending across the South Indian Ocean towards West Australia. These two positive anomalies are coupled by negative anomalies in the Equatorial East Pacific and across South and East Eurasia. An additional negative anomaly is detected in Antarctica. If we consider that a spatial pattern of the sub-lithospheric mantle geoid should mainly manifest sub-lithospheric mantle density heterogeneities, it becomes clear that a long-wavelength pattern in Figure 3 might be affected by large errors especially at locations with a maximum lithospheric thickness. This becomes even more evident when inspecting the sub-lithospheric mantle geoid in Figure 2f, where a maximum lithospheric deepening is clearly manifested in the geoidal geometry. This indicates a possible presence of large errors attributed to lithospheric model uncertainties used in the gravimetric forward modelling. These errors are estimated in the next section.

5. Error analysis

Large errors are expected in the CRUST1.0 and LITHO1.0 models used for a gravimetric forward modelling in this study. As already mentioned in the preceding section, this is particularly evident from pronounced positive sub-lithospheric mantle geoid anomalies (Figure 2f) that are also still partially exhibited in its long-wavelength pattern (Figure 3). These positive anomalies correspond with the largest cratonic lithospheric thickness (the Laurentian Shield in North America, the Amazonian Shield and São Francisco Craton in South America, the West African Craton, the East European and Siberian Cratons and the Baltic Shield in Eurasia, and the West Australian Craton).

As seen in Eq. (47), the lithospheric thickness uncertainties propagate almost linearly to the sub-lithospheric mantle geoid uncertainties. We can, therefore, readily estimate the geoid errors for particular lithospheric thickness uncertainties. According to the LITHO1.0 model, the lithospheric mantle density varies roughly from 3000 to 3450 kg.m^-3. From these density variations we can assume that $\delta {\rho }^{LM}$ is mostly within ±200 kg.m^-3. The boundary between the lowermost lithosphere and the uppermost asthenosphere (i.e., the LAB) is rheological, conventionally taken at the 1300°C isotherm, above which the mantle behaves in a rigid fashion and below which behaves in a ductile fashion [100]. Studies suggest the existence of a compositional or chemical density contrast 0-20 kg.m-3 [101], except probably the cratonic mantle. In addition to a possibly chemical density contrast, a strong thermal density contrast 30-60 kg.m^-3 occurs when the asthenosphere is locally uplifted by rifting. We, therefore, assume that the density variations $\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}$ at the boundary between the lithosphere and the asthenosphere are mostly within ±30 kg.m^-3. For $\delta {\rho }^{LM}+\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}=\pm 230$ kg.m^-3, the errors in the lithospheric thickness estimates of ${\epsilon }_L\approx \pm 10$ km cause errors in the sub-lithospheric mantle geoid roughly ${\epsilon }_{N^{\mathrm{SM}}}\approx \pm 0.1$ km.

The lithospheric thickness models vary significantly. Pasyanos et al. [70], for instance, found large differences between their LITHO1.0 model and that prepared by Conrad and Lithgow-Bertelloni [102]. They also provided explanation that these large differences are due to applying different methods, assumptions, and input data. To get a rough idea about a magnitude of these errors, we compared the lithospheric thickness models SLNAAFSA [103], SL2013sv [104], LITHO1.0 [70], CAM2016 [105,106], and 3D2015-07Sv [107], and used their differences to estimate (sub-lithospheric mantle) geoid errors. We note here that the SLNAAFSA model was prepared by Hoggard et al. [103] by merging regional updates from North America [108], Africa [109], and South America [110] into the global SL2013sv model. The lithospheric thickness models are shown in Figure 4, with their statistical summary given in Table 3. The statistics of differences between individual lithospheric thickness models are summarized in Table 4.

As seen in Table 4, differences between lithospheric thickness models exceed several hundreds of meters. If we consider that average uncertainties in the lithospheric thickness ${\epsilon }_L$ are about ±40 km (the average of the Root-Mean-Square (RMS) of differences in Table 4), the sub-lithospheric mantle geoid errors ${\epsilon }_{N^{\mathrm{SM}}}$ could reach or even exceed ±0.4 km.

Tomographic images of the lithosphere have a limited accuracy and resolution. Moreover, the conversion between seismic wave velocities and rock densities is not unique—e.g., [58,111]. Hager and Richards [112], for instance, mentioned that a lithospheric mantle composition and seismic anisotropy is more complex than a deep mantle structure, making it harder to convert seismic anomalies to density anomalies or buoyancy (cf. also Simmons et al., [113]). We, therefore, expect also large (sub-lithospheric mantle) geoid modelling errors due to lithospheric mantle density uncertainties.

According to Eq. (51), the sub-lithospheric mantle geoid error ${\epsilon }_{N^{\mathrm{SM}}}$ depends on the accuracy of anomalous lithospheric mantle density $\delta {\rho }^{LM}$, and additional errors ${\epsilon }_{\Delta {\rho }^{\mathrm{L}/\mathrm{A}}}$ and ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}$ in adopted values of the LAB density contrast $\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}$ and the Moho density contrast $\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}$. Large density variations occur within the lithospheric mantle, although in comparison to the crust, fewer constraints exist on the composition of the lithospheric mantle. Petrological models for the origin of basalt are used to compile models of the oceanic lithospheric mantle. An undepleted peridotitic source that partially melts to produce the mid-ocean ridge basalt, leaves behind a depleted residue (a harzburgite), the thickness of which depends on the degree of a partial melt. As the oceanic lithosphere cools, the undepleted mantle becomes part of the lithospheric column and its thickness increases with the ocean-floor age. The composition of the continental lithospheric mantle depends on a tectonic province and its age. Jordan [114] proposed that the lithospheric mantle might be cold and buoyant in the thickest cratonic portions of very fast seismic velocity (see also Gung et al., [115]). These conditions require the presence of a depleted layer overlying an undepleted upper mantle layer. Based on the LITHO1.0 lithospheric mantle density variations (roughly within ±200 kg.m^-3) and the previously summarized studies, we can anticipate that lithospheric mantle density uncertainties ${\epsilon }_{\delta {\rho }^{LM}}$ are mostly within ±100 kg.m^-3. We further assume that uncertainties in the average LAB density contrast ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}$ are roughly ±15 kg.m^-3.

Results from seismic and gravimetric studies revealed that the Moho density contrast varies substantially. Artemieva [116] demonstrated that the in situ density contrast between the crystalline crust and the lithospheric upper mantle differs substantially between the Cratonic and the Phanerozoic crust within the European part of the Eurasian tectonic plate, approximately 500 kg.m^-3 for the crust beneath Western Europe, between 300 and 350 kg.m^-3 for most of the Precambrian crust, and by as little as 150-250 kg.m^-3 for some parts of the Baltic and the Ukrainian Shields. Sjöberg and Bagherbandi [117] provided the average values 678±78 and 334±108 kg.m^-3 of the Moho density contrast for the continental and oceanic areas respectively. According to these studies, we speculate that the errors of Moho density contrast ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}$ could be roughly ±50 kg.m^-3. Adopting the above error estimates of ${\epsilon }_{\delta {\rho }^{LM}}\approx \pm 100$ kg.m^-3, ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}\approx \pm 15$ kg.m^-3, and ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}\approx \pm 50$ kg.m^-3 and considering the lithospheric thickness $L$ ~100 km and the Moho depth $M$ of ~30 km, the geoid modelling error ${\epsilon }_{N^{\mathrm{SM}}}$could, according to Eq. (51), reach or even exceed ±0.5 km.

To confirm the above theoretical estimates, we rearranged the error propagation in Eq. (51) into the form that directly incorporates density uncertainties within the uppermost asthenosphere and the lowermost crust. Inserting for $\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}={\mathsf{\rho }}^{\mathrm{LM}}-{\mathsf{\rho }}^{\mathrm{c}}$, $\delta {\rho }^{LM}={\mathsf{\rho }}^{\mathrm{LM}}-{\rho }^{LM}$ and $\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}={\mathsf{\rho }}^{\mathrm{A}}-{\mathsf{\rho }}^{\mathrm{LM}}$ in Eq. (51), we get

$${\epsilon }_{N^{\mathrm{SM}}}\approx -\frac{3}{4\pi {\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}\left[\left({\epsilon }_{{\mathsf{\rho }}^{\mathrm{LM}}}-{\epsilon }_{{\mathsf{\rho }}^{\mathrm{c}}}\right)M+\left({\epsilon }_{{\mathsf{\rho }}^{\mathrm{LM}}}-{\epsilon }_{{\rho }^{LM}}\right)\left(L-M\right)+\left({\epsilon }_{{\mathsf{\rho }}^{\mathrm{A}}}-{\epsilon }_{{\mathsf{\rho }}^{\mathrm{LM}}}\right)L\right]$$

$$\approx -\frac{3}{4\pi {\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}\left[{\epsilon }_{{\mathsf{\rho }}^{\mathrm{LM}}}M-{\epsilon }_{{\mathsf{\rho }}^{\mathrm{c}}}M+{\epsilon }_{{\mathsf{\rho }}^{\mathrm{LM}}}L-{\epsilon }_{{\mathsf{\rho }}^{\mathrm{LM}}}M-{\epsilon }_{{\rho }^{LM}}\left(L-M\right)+{\epsilon }_{{\mathsf{\rho }}^{\mathrm{A}}}L-{\epsilon }_{{\mathsf{\rho }}^{\mathrm{LM}}}L\right]$$

$$\approx \frac{3}{4\pi {\overline{\mathsf{\rho }}}^{\mathrm{Earth}}}\left[{\epsilon }_{{\mathsf{\rho }}^{\mathrm{c}}}M+{\epsilon }_{{\rho }^{LM}}\left(L-M\right)+{\epsilon }_{{\mathsf{\rho }}^{\mathrm{A}}}L\right],$$

(52)

where ${\epsilon }_{{\mathsf{\rho }}^{\mathrm{c}}}$ and ${\epsilon }_{{\mathsf{\rho }}^{\mathrm{A}}}$ are errors in the average crustal ${\mathsf{\rho }}^{\mathrm{c}}$and asthenospheric ${\mathsf{\rho }}^{\mathrm{A}}$ densities respectively. As seen in Eq. (52), the modelling error ${\epsilon }_{N^{\mathrm{SM}}}$ depends on the accuracy of lithospheric mantle density distribution as well as adopted average densities of the crust and the asthenosphere, while errors in adopted value of the average lithospheric mantle density ${\mathsf{\rho }}^{\mathrm{LM}}$ cancel.

Christensen and Mooney [111] reported the average value 2835 kg.m-3 for the continental crustal density. Tenzer et al. [50] provided a very similar value 2790 kg.m-3. The oceanic crust, composed primarily of mafic rocks, is typically heavier than the continental crust—e.g., [118]. Tenzer et al. [50] estimated that the average density of the oceanic crust (without marine sediments) is 2860 kg.m-3. Carlson and Raskin [119] provided the estimate of 2890 kg.m-3. Tenzer and Gladkikh [92] confirmed a similar value (when disregarding marine sediments) of 2900 kg.m-3 based on the analysis of global samples of marine bedrock densities. If we focus only on the lowermost crustal density uncertainties, the density according to the CRUST1.0 varies from 2850 kg.m-3 under Himalayas to 3050 kg.m-3 under oceans with a standard deviation ~70 kg.m-3. Christensen and Mooney [111] provided the uncertainty of ±50 kg.m-3 for the crystalline crust that could be reduced to ~±30 kg.m-3 when considering only the lowermost crust. The error ${\epsilon }_{{\mathsf{\rho }}^{\mathrm{c}}}$ of ±50 kg.m-3 could then be a realistic expectation for an average density within the lowermost crust.

Density variations within the asthenosphere are probably controlled by a mantle flow pattern. Large density variations attributed to mantle plumes, for instance, were detected under Iceland, Hawaii, Azores, and elsewhere—e.g., [120]. Seismic tomography results indicate also relatively large density heterogeneities and density discontinuities in the mantle transition zone at depths approximately between 410 to 660 km—e.g., [121,122]. Based on petrological evidences, Griffin et al. [123] reported density variations 3300-3400 kg.m^-3 in the asthenosphere. From these estimates, we could assume that the error ${\epsilon }_{{\mathsf{\rho }}^{\mathrm{A}}}$ in the average asthenosphere density is roughly ±30 kg.m^-3. Assigning these density uncertainties to Eq. (52), the geoid modelling errors ${\epsilon }_{N^{\mathrm{SM}}}$ within ±0.6 km closely corresponds with the estimated value of ±0.5 km according to Eq. (51).

To check the correctness of our initial assumption that the lithospheric mantle density uncertainties ${\epsilon }_{\delta {\rho }^{LM}}$ are mostly within ±100 kg.m^-3, we used differences between the lithospheric thickness models (in Table 4) to reversely estimate errors of lithospheric mantle density. The computation was done approximately by using a simple relation for the gravitational potential $V$ of an infinite planar plate that is computed from the density $\rho$ and thickness $s$ values as follows

$$V\approx 2\pi GR\rho s$$

(53)

The error propagation model for Eq. (53) reads

$${\epsilon }_V\approx 2\pi GR\left(s{\epsilon }_{\rho }+\rho {\epsilon }_s\right),$$

(54)

where ${\epsilon }_{\rho }$ is a density uncertainty, and ${\epsilon }_s$ is a thickness uncertainty.

Since we here inspect the propagation only between errors ${\epsilon }_{\delta {\rho }^{LM}}$ and ${\epsilon }_{N^{\mathrm{SM}}}$ (or equivalently ${\epsilon }_V$), the error ${\epsilon }_s$ in Eq. (54) is disregarded. We then write

$${\epsilon }_V\approx 2\pi GRs{\epsilon }_{\rho }$$

(55)

Defining the potential errors ${\epsilon }_V$ in Eq. (55) in terms of geoid uncertainties ${\epsilon }_N$ , we arrive at

$${\epsilon }_N=\frac{{\epsilon }_V}{{\gamma }_o}{\epsilon }_V\approx \frac{2\pi }{{\gamma }_o}GRs{\epsilon }_{\rho }\approx 2\pi \frac{GR}{\mathrm{GM}/R^2}s{\epsilon }_{\rho }\approx 2\pi \frac{R^3}{\mathrm{M}}s{\epsilon }_{\rho },$$

(56)

where M = 5.9722 x 10²⁴ kg is the Earth’s total mass.

According to Table 4, the RMS of differences between the lithospheric thickness models is on average roughly ±0.4 km. This RMS corresponds to the density error ${\epsilon }_{\rho }$ of ~ ±160 kg.m^-3. As seen, this error estimate agrees quite closely with the error analysis according to Eq. (51), where we assumed density errors within a similar range (i.e., ${\epsilon }_{\delta {\rho }^{LM}}\approx \pm 100$ kg.m^-3, ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{L}/\mathrm{A}}}\approx \pm 15$ kg.m^-3, and ${\epsilon }_{\Delta {\mathsf{\rho }}^{\mathrm{c}/\mathrm{m}}}\approx \pm 50$ kg.m^-3).

As confirmed in this study, the sub-lithospheric mantle geoid modelling is affected mainly by errors in lithospheric thickness and lithospheric mantle density. Since we disregarded density variations within the asthenosphere, additional errors in our result are possibly attributed to temperature anomalies and small-scale convection patterns within the asthenosphere. It is also clear from Eq. (52) that the geoid modelling errors due to lithospheric mantle density uncertainties increase almost linearly with the increasing lithospheric thickness. Consequently, both uncertainties likely propagate into maximum geoid errors at locations characterized by the largest continental lithospheric deepening under cratonic formations that are manifested by large positive anomalies of the sub-lithospheric mantle geoid model (Figure 2e). At these locations, errors could according to our estimates reach (or even exceed) ±0.5 km.

These theoretical findings indicate that the sub-lithospheric mantle geoid presented in Figure 2f (or Figure 3) might be very inaccurate, especially at locations with a maximum lithospheric deepening where these errors reach maxima. In reality, modifications of the geoidal geometry by sub-lithospheric mantle density heterogeneities are probably much smaller than those seen in Figure 2f. This is due to the fact that relative lateral density variations within the mantle are up to only a few percent—e.g., [124,125,126,127,128,129].

6. Summary and concluding remarks

We have analyzed errors of computing the sub-lithospheric mantle geoid based on applying a gravimetric forward modelling and by using lithospheric density and thickness models as input data. The principle of this method is to model and subsequently remove the gravitational contributions of topography and lithospheric density and geometry variations in order to enhance a gravitational signature of the sub-lithospheric mantle density structure in the geoidal geometry. The expressions in the spectral domain, that utilize methods for a spherical harmonic analysis and synthesis of gravitational and lithospheric density structure models, were then used to establish the error propagation between geoid modelling errors and lithospheric mantle and lithospheric thickness uncertainties.

According to our analysis, errors in modelling of sub-lithospheric mantle geoid due to lithospheric thickness and lithospheric mantle density uncertainties could reach or even exceed ±0.5 km, particularly at locations with a maximum lithospheric thickness under cratons where both errors could magnify significantly. Theoretical error analysis and numerical results, therefore, suggest that actual modifications of the geoidal geometry by lateral mantle density variations below the lithosphere-asthenosphere boundary are much smaller than those computed and presented in Figure 2f (or in Figure 3).

This finding is indicative for a possible improvement of the Earth’s synthetic density and gravitational models based on a density calibration within individual lithospheric and underlying deeper mantle layers so that a sum of volumetric mass densities within the Earth’s synthetic model provides the value that closely agrees with the Earth’s total mass (excluding the atmosphere). Equivalently, the Earth’s synthetic density model should generate the gravitational field that closely agree with the Earth’s gravitational field. This issue will be addressed in the forthcoming study, where we will compile a synthetic model that closely fulfills both requirements.