Sensing the global CRUST1.0 Moho by gravitational curvatures of crustal mass anomalies

ABSTRACT

The gravitational curvatures (GC, the third-order derivatives of the gravitational potential) in gravity field modeling are gaining increased interest in geosciences. The crustal effects of the GC and Moho variation sensed by the GC are not fully evaluated for the current study. In this contribution, the effects of the GC induced by topographic and anomalous crustal masses are investigated based on ETOPO1 and CRUST1.0 models using the tesseroids. By adopting the gravitational stripping correction, the residual GC sensed by the CRUST1.0 Moho depths are presented globally to examine whether the GC can sense crustal mass anomalies at the satellite altitude of 260 km . The spatial analysis using the Pearson correlations coefficient (PCC) between the residual GC and the CRUST1.0 Moho depths is performed. Among the 10 residual GC functionals, the PCC value of the residual radial-radial-radial component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ is largest with 0.31, where this value is highly dependent on the spectral content removed from the EGM2008, e.g. signals assumed to relate to deeper mass anomalies. Numerical experiments show that with the increased order of the derivatives up to third-order, the fineness level of different global Moho sensed crustal mass anomalies increases correspondingly. Taking the Tibetan plateau for example, the values of the $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ can better reveal the detailed features of the Tibetan plateau’s Moho depth than these of the lower-order residual radial functionals (i.e. disturbing potential $\mathit{\delta }{T}^{\mathrm{res}}$, disturbing radial gravity vector $\mathit{\delta }{T}_z^{\mathrm{res}}$, and disturbing radial-radial gravity gradient tensor $\mathit{\delta }{T}_{\mathit{zz}}^{\mathrm{res}}$), especially for the Qaidam, Sichuan, Tarim, and Turpan basins. Numerical results over the Himalayan region demonstrate that the GC component $\mathit{\delta }{T}_z^{\mathrm{res}}$ has some potential in geophysical inversion. These residual GC functionals would help to get a better knowledge of the internal structures of the Earth and other planetary objects.

KEYWORDS:

• Gravity field
• gravitational curvatures
• tesseroid
• gravitational stripping correction
• Moho

1. Introduction

In the field of gravity forward modeling, the gravitational stripping correction is one of the important techniques on studying the internal structure of the Earth (Tenzer, Hamayun, and Vajda 2009a, 2009b; Novák and Tenzer 2013; Tenzer et al. 2015; Chen 2017; Chen, Tenzer, and Li 2018). The definition of the gravitational stripping correction can be described to perform the correction on the observed Earth gravity field based on the Earth gravity model by stripping the different gravitational signals of the internal mass anomalies of the Earth. The gravitational stripping correction was the process of starting from the Earth gravity model minus the gravitational signals assumed to be known mass distributions. The digital elevation models and crustal models were used to model the gravitational signals of crustal mass anomalies. After the gravitational stripping correction, the residual gravitational signals were obtained to reveal the internal mass variations of the gravitational signals.

Various studies were performed on the topographic effects, crustal effects, and the residual gravitational components’ variations with the Mohorovičić discontinuity (Moho) depths after the gravitational stripping correction. Tenzer, Hamayun, and Vajda (2009a) analyzed the correlation of the residual gravity disturbances, gravity anomalies and disturbing potentials with the topography, bathymetry from the ETOPO5 model (e.g. Edwards (1989) with a grid size of $5^{{ }^{\circ }}\times 5^{{ }^{\circ }}$) and the Moho depths from the CRUST2.0 model (e.g. Bassin, Laske, and Masters (2000) with a grid size of $2^{{ }^{\circ }}\times 2^{{ }^{\circ }}$). Tenzer, Hamayun, and Vajda (2009b) evaluated the global topographic and crustal effects of gravity disturbances with the CRUST2.0 model. Novák (2010) evaluated the gravitational potential (GP) values of the solid topography, seawater and atmosphere in the spectral domain, where the forms of the formulae were represented in the spherical harmonic coefficients. Tenzer et al. (2012) investigated the spatial and spectral properties of the four residual gravity disturbances (i.e. consolidated crust-stripped, complete crust-corrected, homogeneous crust, and crust-stripped relative-to-mantle) from the CRUST2.0 model in the spectral domain. Novák et al. (2013) calculated the crustal effects of the disturbing gravitational gradients from the topography, continental ice, and layers (i.e. soft and hard sediments, upper, middle and lower crust) of the CRUST2.0 model in the spectral domain. Novák and Tenzer (2013) performed the residual disturbing gravitational gradients (e.g. $\mathit{\delta }{\mathrm{\Gamma }}_{\mathit{xx}}$, $\mathit{\delta }{\mathrm{\Gamma }}_{\mathit{yy}}$, and $\mathit{\delta }{\mathrm{\Gamma }}_{\mathit{zz}}$) variations with the CRUST2.0 Moho depths. Tenzer et al. (2015) evaluated the residual gravity disturbances variations ($\mathit{\delta g}$) with the CRUST1.0 Moho depths (e.g. Laske et al. (2013) with a grid size of $1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$). Uieda and Barbosa (2017) calculated crustal effects of gravity disturbance from the ETOPO1 (e.g. Amante and Eakins (2009) with a grid size of $1^{\prime }\times 1^{\prime }$) and CRUST1.0 model for South America using tesseroids, then applied the nonlinear gravity inversion to obtain the Moho depths of the South America. Chen (2017) performed the gravity forward modeling with the residual vertical gravity gradient tensor (GGT, i.e. the second-order derivatives of the GP) using tesseroids and obtained the Moho depth by gravimetric inversion with the residual gravity-gradient data in the Tibet region. Chen, Tenzer, and Li (2018) calculated the residual GP by subtracting crustal effects of the GP from the GP model values and performed the gravimetric inversion to reveal the Moho depths of the Tibet area. Deng and Shen (2019) investigated the topographic effects of the gravitational curvatures (GC, i.e. the third-order derivatives of the GP) with the ETOPO1 over China. Novák et al. (2019) presented that third-order gradients were useful for the evaluation of the shallow structures and subsurface with the radial-radial-radial third-order gradients of the upper sediments from the CRUST1.0 model. Regarding the recent research advances of the GC, Romeshkani, Sharifi, and Tsoulis (2020) investigated the joint recovery of gravity anomalies with the GGT and GC including the usage of the topographic/isostatic gravity model. Romeshkani, Sharifi, and Tsoulis (2021) performed the estimation of the GC using the GGT from the Gravity field and steady-state Ocean Circulation Explorer mission (GOCE) data.

From the above literature review, current studies of the topographic effects, crustal effects, and the sensed Moho variations mainly focused on the lower-order gravitational functionals, i.e. GP (Novák 2010; Chen, Tenzer, and Li 2018), gravity disturbances (Tenzer, Hamayun, and Vajda 2009a, 2009b; Tenzer et al. 2012, 2015; Uieda and Barbosa 2017), gravity anomalies (Tenzer, Hamayun, and Vajda 2009a), and GGT (Novák et al. 2013; Novák and Tenzer 2013; Chen 2017), whereas the crustal effects of the GC and Moho variation sensed by the GC are not fully evaluated (Deng and Shen 2019; Novák et al. 2019). The Moho is a smooth surface, and the crustal rigidity compensate loads, where even the GGT components have been applied to sense the Moho variations (Novák and Tenzer 2013; Chen 2017) and it means that the second-order derivatives contain high frequencies as well. Thus, it is essential to investigate the crustal effects of the GC and find out whether the GC could sense the Moho variations, which will help to clarify whether the GC can be used in Moho’s related research, especially the inversion of the internal structure of the Earth. This paper intends to fill these research gaps and deals with the crustal effects and Moho variations sensed by the GC functionals compared to those of the lower-order gravitational functionals (i.e. GP, gravity vector (GV, i.e. the first-order derivatives of the GP), and GGT).

For the gravity modeling techniques, there are spectral modelings to represent the Earth gravity model (i.e. in the form of the spherical harmonic coefficients) and spatial gravity forward modeling in this contribution. Based on the digital elevation and crustal models, there are various mass bodies to model the topographic and crustal effects for gravity forward modeling in the spatial domain, for instance, linear vertical mass, layer mass, point mass, prism, polyhedron, tesseroid, spherical triangular tessellation, and arbitrary shape mass (Nagy, Papp, and Benedek 2000; Heck and Seitz 2007; Wild-Pfeiffer 2008; Tsoulis 2012; Zhang, Mooney, and Chen 2018; Wu 2018a; Tsoulis and Gavriilidou 2021). Among these mass bodies, the tesseroid can be considered the elementary mass body, which is recommended to evaluate the topographic effects globally at the satellite height (Uieda, Valéria, and Braitenberg 2016; Deng and Shen 2019). Because its computation efficiency is globally good enough (Heck and Seitz 2007; Wild-Pfeiffer 2008; Grombein, Seitz, and Heck 2013; Deng et al. 2016), and the definition of the tesseroid grid elements in latitude and longitude coincides with these of the released digital elevation and crustal models. Even though there occur large approximation errors at the very near area not only for the GP, GV and GGT functionals (Heck and Seitz 2007; Uieda, Valéria, and Braitenberg 2016; Shen and Deng 2016), but also for the GC functionals (Deng and Shen 2018a, 2018b, 2019), its relative approximation errors (about 0.1%) are acceptable for the GP, GV, GGT and GC functionals at the satellite height (Uieda, Valéria, and Braitenberg 2016; Deng and Shen 2019). Based on these advantages, the tesseroids are applied to model the topographic and crustal effects of the GC in this paper, which is the continuous work of Deng and Shen (2019).

In this contribution, the GC effects from topographic and crustal mass anomalies are investigated with the ETOPO1 and CRUST1.0 models using the tesseroids, where the crustal effects in the definition of this paper include the topographic effects of the disturbing gravitational components. After performing the gravitational stripping correction, the spatial correlations between the residual GC functionals and the CRUST1.0 Moho depths are calculated using the Pearson correlation coefficients. We investigate the un-modeled effects of the degree/order (d/o) value on the Pearson correlation coefficients, where the un-modeled effects are the deeper mass anomalies in the mantle and core, and these have their spectral energy at a very low d/o value. Moreover, the Moho variations sensed by the GC at the Himalayan region are performed in comparison with the low-order gravitational functionals.

The paper is organized as follows. Section 2 presents the theoretical methodologies in which Section 2.1 presents the spectral modeling for an Earth gravity model of the GC. In Section 2.2, the spatial modeling for the GC effects of the tesseroid from crustal mass distributions is presented. Section 2.3 provides the formula of the residual GC functionals with gravitational stripping correction. In Section 2.4, the formula of Pearson correlation coefficient for spatial correlation analysis is reviewed. The numerical experiments are performed in Section 3. Section 3.1 provides the setup of the numerical experiments. In Section 3.2, the crustal effects of the GC functionals are calculated with the ETOPO1 and CRUST1.0 models. Section 3.3 presents the residual GC functionals after the gravitational stripping correction. The Laplace parameters of the GC are investigated in Section 3.4. In Section 3.5, the spatial analysis with Pearson correlation coefficient between the residual GC functionals and the CRUST1.0 Moho depths is performed. The comparison between the residual GC functionals and the lower-order gravitational functionals sensed by the Moho depth in the Himalayan region is revealed in Section 3.6. The conclusions of this paper are summarized, and the viewpoints on further research work are recommended in Section 4.

2. Methodology

2.1. Spectral modeling for an Earth gravity model of the gravitational curvatures

In the spectral domain, the radial-radial-radial component of the disturbing GC (denoted as $T_{\mathit{zzz}}^{\mathrm{EGM}}$) can be expressed in spherical harmonic series in the local north-oriented reference frame as (Hamáčková et al. 2016):

(1) $\begin{array}{rl}& T_{\mathit{zzz}}^{\mathrm{EGM}}(\mathit{r},\mathit{\phi },\mathit{\lambda })=-\frac{\mathit{GM}}{R^4}\sum _{\mathit{n}=0}^N{(\frac{R}{r})}^{\mathit{n}+4}(\mathit{n}+1)(\mathit{n}+2)(\mathit{n}+3)\\ & \sum _{\mathit{m}=0}^n\left({\bar{C}}_{\mathit{n},\mathit{m}}\cos \mathrm{m}\mathit{\lambda }+{\bar{S}}_{\mathit{n},\mathit{m}}\mathit{sin}\mathit{m\lambda }\right){\bar{P}}_{\mathit{n},\mathit{m}}(sin\mathit{\phi })\end{array}$(1)

where $\mathit{GM}$ is the geocentric gravitational constant. $\mathit{R}$ is the mean Earth radius, where the basic theory of the reference sphere can be referred to Kuhn (2000) and Tsoulis (1999, 2004). $\mathit{r}$, $\mathit{\phi }$, and $\mathit{\lambda }$ are the spherical radius, latitude and longitude of the computation point (or field point). ${\bar{P}}_{\mathit{n},\mathit{m}}(sin\mathit{\phi })$ is the first kind fully normalized and associated Legendre function with degree $\mathit{N}$ and order $\mathit{m}$ ($\mathit{N}$ is the maximum degree). ${\bar{C}}_{\mathit{n},\mathit{m}}$ and ${\bar{S}}_{\mathit{n},\mathit{m}}$ are the fully normalized spherical harmonic coefficients in Equation (1)(1) $\begin{array}{rl}& T_{\mathit{zzz}}^{\mathrm{EGM}}(\mathit{r},\mathit{\phi },\mathit{\lambda })=-\frac{\mathit{GM}}{R^4}\sum _{\mathit{n}=0}^N{(\frac{R}{r})}^{\mathit{n}+4}(\mathit{n}+1)(\mathit{n}+2)(\mathit{n}+3)\\ & \sum _{\mathit{m}=0}^n\left({\bar{C}}_{\mathit{n},\mathit{m}}\cos \mathrm{m}\mathit{\lambda }+{\bar{S}}_{\mathit{n},\mathit{m}}\mathit{sin}\mathit{m\lambda }\right){\bar{P}}_{\mathit{n},\mathit{m}}(sin\mathit{\phi })\end{array}$(1) of Hamáčková et al. (2016). The definition of the local north-oriented reference frame is a system with its origin at the field point, and its x-, y-, and z-axis point to the north, east, and up along the radial direction, respectively. The other disturbing GC functionals (e.g. $T_{\mathit{xxx}}^{\mathrm{EGM}}$, $T_{\mathit{xxy}}^{\mathrm{EGM}}$, $T_{\mathit{xxz}}^{\mathrm{EGM}}$, $T_{\mathit{xyz}}^{\mathrm{EGM}}$, $T_{\mathit{yyx}}^{\mathrm{EGM}}$, $T_{\mathit{yyy}}^{\mathrm{EGM}}$, $T_{\mathit{yyz}}^{\mathrm{EGM}}$, $T_{\mathit{zzx}}^{\mathrm{EGM}}$, and $T_{\mathit{zzy}}^{\mathrm{EGM}}$) in the local north-oriented reference frame have been derived in Equations (22) – (30) of Hamáčková et al. (2016).

In this study, an Earth gravity model with spectral modeling is applied in the next numerical experiments as the initial values in the gravitational stripping correction. Different Earth gravity models can be chosen at the global gravity field models from the website (http://icgem.gfz-potsdam.de/tom_longtime) of the international center for global Earth models (Ince et al. 2019). Note that the spectral modeling technique is not the focus but only a way to represent the Earth gravity model.

2.2 Spatial modeling for the gravitational curvatures of the tesseroid from crustal mass distributions

In the spatial domain, the radial-radial-radial component of the GC (denoted as $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{tess}}$) of a tesseroid (see Figure 1 of Deng and Shen (2019)) using the Cartesian integral kernels can be obtained in the local north-oriented reference frame as (Deng and Shen 2018b, 2019) :

Figure 1. Crustal effects of the $\mathit{\delta }{T}_{\mathit{zzz}}$ values based on the ETOPO1 and CRUST1.0 by a grid size of $1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$ at the satellite height of 260 km: (a) topography/bathymetry effects, (b) ice effects, (c) upper sediment effects, (d) middle sediment effects, (e) lower sediment effects, (f) upper consolidated crust effects, (g) middle consolidated crust effects, and (h) lower consolidated crust effects.

(2) $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{tess}}\left(\mathit{r},\mathit{\phi },\mathit{\lambda }\right)=\mathit{G}\mathrm{\Delta }\mathit{\rho }\underset{{\lambda }_1}{\overset{{\lambda }_2}{\int }}\underset{{\phi }_1}{\overset{{\phi }_2}{\int }}\underset{r_1}{\overset{r_2}{\int }}{{\mathit{r}}^{\prime }}^{2}cos\mathit{\phi }{\mathit{\hspace{0.17em}}}_{\mathit{\hspace{0.17em}}}^{\prime }\left(\frac{15{\mathrm{\Delta }}_z^3}{{\mathrm{\ell }}^7}-\frac{9{\mathrm{\Delta }}_z}{{\mathrm{\ell }}^5}\right)\mathit{dr}{\mathit{\hspace{0.17em}}}_{\mathit{\hspace{0.17em}}}^{\prime }\mathit{d\phi }{\mathit{\hspace{0.17em}}}_{\mathit{\hspace{0.17em}}}^{\prime }\mathit{d}{\lambda }^{\prime }$(2)

(3) ${\mathrm{\Delta }}_z=r^{\prime }\mathit{cos\psi }-\mathit{r}$(3)

(4) $\mathrm{\ell }=\sqrt{{r^{\prime }}^2+r^2-2r^{\prime }\mathit{r}cos\mathit{\psi }}$(4)

(5) $cos\mathit{\psi }=sin\mathit{\phi }sin{\phi }^{\prime }+cos\mathit{\phi }cos{\phi }^{\prime }cos\left({\lambda }^{\prime }-\mathit{\lambda }\right)$(5)

where $\mathit{G}$ is the Newtonian gravitational constant. $\mathrm{\Delta \rho }$ is the density contrast with respect to the reference crust density ${\rho }^{\mathit{crust}}$. $r^{\prime }$, ${\phi }^{\prime }$, and ${\lambda }^{\prime }$ are the spherical radius, latitude and longitude of the integration point. In addition, the formulae of other 9 disturbing GC functionals (e.g. $\mathit{\delta }{T}_{\mathit{xxx}}^{\mathrm{tess}}$, $\mathit{\delta }{T}_{\mathit{xxy}}^{\mathrm{tess}}$, $\mathit{\delta }{T}_{\mathit{xxz}}^{\mathrm{tess}}$, $\mathit{\delta }{T}_{\mathit{xyz}}^{\mathrm{tess}}$, $\mathit{\delta }{T}_{\mathit{yyx}}^{\mathrm{tess}}$, $\mathit{\delta }{T}_{\mathit{yyy}}^{\mathrm{tess}}$, $\mathit{\delta }{T}_{\mathit{yyz}}^{\mathrm{tess}}$, $\mathit{\delta }{T}_{\mathit{zzx}}^{\mathrm{tess}}$, and $\mathit{\delta }{T}_{\mathit{zzy}}^{\mathrm{tess}}$) in the local north-oriented reference frame can be referred to in Equation (1)(1) $\begin{array}{rl}& T_{\mathit{zzz}}^{\mathrm{EGM}}(\mathit{r},\mathit{\phi },\mathit{\lambda })=-\frac{\mathit{GM}}{R^4}\sum _{\mathit{n}=0}^N{(\frac{R}{r})}^{\mathit{n}+4}(\mathit{n}+1)(\mathit{n}+2)(\mathit{n}+3)\\ & \sum _{\mathit{m}=0}^n\left({\bar{C}}_{\mathit{n},\mathit{m}}\cos \mathrm{m}\mathit{\lambda }+{\bar{S}}_{\mathit{n},\mathit{m}}\mathit{sin}\mathit{m\lambda }\right){\bar{P}}_{\mathit{n},\mathit{m}}(sin\mathit{\phi })\end{array}$(1) and (A.11) – (A.19) of Deng and Shen (2019).

The total GC effects of the respective mass distributions (i.e. topographic and crustal models) are obtained as:

(6) $\mathit{\delta }{T}_{\mathit{ijk}}=\sum \mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{tess}}$(6)

where $\mathit{i},\mathit{j},\mathit{k}=\mathit{x},\mathit{y},\mathit{z}$ and the total GC effects ($\mathit{\delta }{T}_{\mathit{ijk}}$) are the sum of the individual effects calculated from the tesseroid mass bodies, and $\sum$ means the sum of the GC effects of the tesseroids forming the whole crustal layers (e.g. topography/bathymetry, ice, sediments, and consolidated crusts).

Using the GC formulae of the tesseroid, the latitude and longitude of the individual tesseroid can correspond to those of the grid size of the ETOPO1 and CRUST1 models. Moreover, the computation efficiency and relative approximation errors are sufficient in the numerical operation, which has been confirmed for the topographic effects of the GC at the satellite height of 260 km (Deng and Shen 2019). The difference between this paper and Deng and Shen (2019) lies in that the range of applications has expanded from topographic effects (i.e. topography/bathymetry) to crustal effects (i.e. topography/bathymetry, ice, sediments, and consolidated crusts) of the GC functionals.

2.3. Residual gravitational curvatures sensing the Moho variations after the gravitational stripping correction

The formula of the residual GC functionals sensing the Moho variations after the gravitational stripping correction (denoted as $\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{res}}$) can be expressed as:

(7) $\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{res}}=T_{\mathit{ijk}}^{\mathrm{EGM}}-\{T_{\mathit{ijk}}^{\mathrm{EGM}}{\}}^{\mathit{un}-\mathit{modeled}}-\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{Topo}/\mathrm{bath}}-\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{Ice}}-\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{Sed}}-\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{Con}-\mathrm{crust}}$(7)

where $\mathit{i},\mathit{j},\mathit{k}=\mathit{x},\mathit{y},\mathit{z}$. The superscript EGM represents the formula of the GC functionals with spherical harmonic series in the spectral domain mentioned in Section 2.1. The superscripts Topo/bath, Ice, Sed, and Con-crust mean the different crustal effects using tesseroids in the spatial domain in Section 2.2, which are listed as topography/bathymetry, ice, (upper, middle and lower) sediments, and (upper, middle and lower) consolidated crusts in the ETOPO1 and CRUST1.0 models.

In this study, it is assumed that the gravitational signals of deeper mass anomalies are solely caused by the mantle/core mass anomalies, please refer to Bowin (1983), where the depth of mass anomalies (or point masses) was related to the spectral signal for the gravity anomaly. After neglecting the crustal uncertainties, the superscript un−modeled means the deep Earth layers’ effects from the mantle to the inner core of the Earth.

Equation (7)(7) $\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{res}}=T_{\mathit{ijk}}^{\mathrm{EGM}}-\{T_{\mathit{ijk}}^{\mathrm{EGM}}{\}}^{\mathit{un}-\mathit{modeled}}-\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{Topo}/\mathrm{bath}}-\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{Ice}}-\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{Sed}}-\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{Con}-\mathrm{crust}}$(7) provides the main principle of gravitational stripping correction as: (1) based on the Earth gravity models from the website (http://icgem.gfz-potsdam.de/tom_longtime), the observed (or true) gravity field is obtained; (2) the gravitational influence of the mantle and core below the Moho surface is subtracted from the model values obtained in the first step, where this step can generally be achieved by the spectral filtering technique (Eckhardt 1983; Bowin 1983, 2000; Bagherbandi and Sjöberg 2012; Chen 2017; Chen, Tenzer, and Li 2018), i.e. subtracting the first d/o value 17 spherical harmonic coefficients from the model values; (3) the gravitational effects of the different crustal layers using tesseroids are calculated in the spatial domain, then these effects are subtracted from the values in the second step. After the above three steps, the remaining signal with Moho variations is residual to the observed (or true) gravity field.

2.4. Spatial correlation analysis using the Pearson correlation coefficient

The Pearson correlation coefficient (PCC), denoted as $\mathit{\gamma }$, can describe the linear correlation between two vectors $\vec{X}=\{x_1,\dots ,x_n\}$ and $\vec{Y}=\{y_1,\dots ,y_n\}$. The PCC formula can be given as (Pearson 1895):

(8) $\mathit{\gamma }=\frac{{\sum }_{\mathit{i}=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{{\sum }_{\mathit{i}=1}^n{(x_i-\bar{x})}^2}\sqrt{{\sum }_{\mathit{i}=1}^n{(y_i-\bar{y})}^2}}$(8)

where $\mathit{n}$ is the size, $x_i$ and $y_i$ are the individual sample points indexed with $\mathit{i}=1,\dots ,\mathit{n}$. $\bar{x}=\frac{1}{n}{\sum }_{\mathit{i}=1}^n{x}_i$ and $\bar{y}=\frac{1}{n}{\sum }_{\mathit{i}=1}^n{y}_i$ are the sample mean values of the two vectors $\vec{X}$ and $\vec{Y}$. $\mathit{\gamma }\in [-1,1]$. $\mathit{\gamma }$ = 0 means that $\vec{X}$ and $\vec{Y}$ have no linear correlation; $\mathit{\gamma }$ = 1 or −1 mean that $\vec{X}$ and $\vec{Y}$ have total positive or negative linear correlation, respectively.

The PCC was adopted for the spatial correlation analysis between the residual gravity disturbances and the Moho and topography depths in Tenzer, Hamayun, and Vajda (2009a), Tenzer et al. (2012), and Tenzer et al. (2015). Herein, we apply the PCC to numerically analyze the relationship between the residual GC functionals and Moho depths in Section 3.5. Note that there exist some other methods to compare identical spatial fields that were obtained by some different parameters’ selection, for instance, different Moho depths evaluated for different compensation levels, or use some of the standard existing spectral or spatial measures (Tsoulis and Patlakis 2013), which will be investigated in the future research.

3. Numerical investigations

3.1. Setup of the numerical experiments

The numerical values of the parameters are listed in Table A1. The height above the reference sphere ($\mathit{h}$) is the computation height as $\mathit{h}$ = 260 km. The grid size $1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$ is for topography and crustal mass distributions. The calculations of the TEGM are conducted with the codes from Hamáčková et al. (2016). The standard forward column method (Holmes and Featherstone 2002) is adopted for computing ${\bar{P}}_{\mathit{n},\mathit{m}}(sin\mathit{\phi })$ in Equation (1)(1) $\begin{array}{rl}& T_{\mathit{zzz}}^{\mathrm{EGM}}(\mathit{r},\mathit{\phi },\mathit{\lambda })=-\frac{\mathit{GM}}{R^4}\sum _{\mathit{n}=0}^N{(\frac{R}{r})}^{\mathit{n}+4}(\mathit{n}+1)(\mathit{n}+2)(\mathit{n}+3)\\ & \sum _{\mathit{m}=0}^n\left({\bar{C}}_{\mathit{n},\mathit{m}}\cos \mathrm{m}\mathit{\lambda }+{\bar{S}}_{\mathit{n},\mathit{m}}\mathit{sin}\mathit{m\lambda }\right){\bar{P}}_{\mathit{n},\mathit{m}}(sin\mathit{\phi })\end{array}$(1) as well as for all remaining GC functionals as given by Hamáčková et al. (2016). Meanwhile, the EGM2008 model (Pavlis et al. 2012) is chosen for the Earth gravity field as the initial value.

The calculation of the crustal effects for the GC functionals is conducted in Section 3.2. The data from the global ETOPO1 digital terrain model (Amante and Eakins 2009) and CRUST1.0 crustal model (Laske et al. 2013) is used by a grid size of $1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$ for the modeling of the 8 crustal compositions. The density of the topography (${\rho }^{\mathrm{topo}}$) and density-contrast of the bathymetry ($\mathit{\delta }{\rho }^{\mathrm{bath}}$) in the ETOPO1 model are listed in Table A1, and the density-contrasts of other crustal layers in the CRUST1.0 model are obtained by subtracting the reference crust density ${\rho }^{\mathrm{crust}}$ from the corresponding density values of each layer’s element in CRUST1.0, which is the same method as used in Uieda and Barbosa (2017).

In Section 3.3, the gravitational stripping correction of the GC functionals are presented with respect to the EGM2008 model and the results of crustal effects in Section 3.2. After applying all gravitational stripping corrections, it is assumed that the residual signal is only the mass anomalies due to the Moho. Subsequently, the residual signal is compared to the CRUST1.0 Moho depths. After the gravitational stripping correction, the variations of the global residual GC functionals sensed with the CRUST1.0 Moho depths are obtained.

In Section 3.5, we perform the spatial analysis using the PCC values between the residual GC functionals and the CRUST1.0 Moho depths. In addition, the effects of the selected d/o value for the spectral filtering technique (Eckhardt 1983; Bagherbandi and Sjöberg 2012; Chen 2017; Chen, Tenzer, and Li 2018) in the spectral domain on the PCC values are analyzed.

3.2. The GC effects of crustal mass anomalies

In order to examine whether the GC functionals ($\mathit{\delta }{T}_{\mathit{ijk}}$) could be used to identify crustal mass anomalies, the crustal effects of the GC functionals are derived based on the methodology in Section 2.2. Calculating different crustal layers with tesseroids, the adaptive discretization stack-based algorithm by the Gauss-Legendre quadrature method is applied for $\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{tess}}$ in Equation (2)(2) $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{tess}}\left(\mathit{r},\mathit{\phi },\mathit{\lambda }\right)=\mathit{G}\mathrm{\Delta }\mathit{\rho }\underset{{\lambda }_1}{\overset{{\lambda }_2}{\int }}\underset{{\phi }_1}{\overset{{\phi }_2}{\int }}\underset{r_1}{\overset{r_2}{\int }}{{\mathit{r}}^{\prime }}^{2}cos\mathit{\phi }{\mathit{\hspace{0.17em}}}_{\mathit{\hspace{0.17em}}}^{\prime }\left(\frac{15{\mathrm{\Delta }}_z^3}{{\mathrm{\ell }}^7}-\frac{9{\mathrm{\Delta }}_z}{{\mathrm{\ell }}^5}\right)\mathit{dr}{\mathit{\hspace{0.17em}}}_{\mathit{\hspace{0.17em}}}^{\prime }\mathit{d\phi }{\mathit{\hspace{0.17em}}}_{\mathit{\hspace{0.17em}}}^{\prime }\mathit{d}{\lambda }^{\prime }$(2) (Deng and Shen 2019). The computation height $\mathit{h}$ is 260 km (c.f. Table A1) for future application of the GC, which might be in future satellite missions, and the distance-size ratio value D = 10 is selected for the evaluation of the $\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{tess}}$ (Deng and Shen 2019). In this way, at the computation height of 260 km, the relative approximation error is about 0.01% for $\mathit{\delta }{T}_{\mathit{ijk}}^{\mathrm{tess}}$ from Figure 5 of Deng and Shen (2019), which is below the threshold error 0.1%. The precision level is satisfactory for geophysical applications.

The numerical calculations of the GC functionals ($\mathit{\delta }{T}_{\mathit{ijk}}$) using tesseroids for the different 8 crustal effects are performed globally. For brevity, only the crustal effects of the GC component $\mathit{\delta }{T}_{\mathit{zzz}}$ are orderly presented in Figure 1. And the statistical information of the GC component $\mathit{\delta }{T}_{\mathit{zzz}}$ for the 8 crustal effects with minimum (Min), maximum (Max), mean (Mean), and standard deviation (STD) are listed in Table A2.

3.3. The residual gravitational curvatures after the gravitational stripping correction

According to the principle of gravitational stripping correction in Section 2.3, we apply the process of gravitational stripping correction beginning with the EGM2008 model as an example to show the values of the residual GC functionals. The residual GC functionals mean the residual values after the different stripping crustal correction from the EGM2008 values. Because of the algorithmic limitations in the standard forward column method (Holmes and Featherstone 2002; Hamáčková et al. 2016) used for the evaluation of the GC functionals in the spectral domain, the truncated maximum degree of the EGM2008 model is 1800. Other information (e.g. the un-modeled effects and the order of the gravitational stripping correction) can be related to the methodology introduced in Section 2.3.

After performing the gravitational stripping correction, the residual GC functionals sensing with the Moho interface are obtained. For brevity, we only demonstrate all residual GC functionals after the gravitational stripping correction in Figure 2, and the statistical information of the GC functionals with minimum (Min), maximum (Max), mean (Mean), and standard deviation (STD) are listed in A3.

Figure 2. The residual GC functionals after the gravitational stripping correction by a grid size of $1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$ at the satellite height of 260 km: (a) $\mathit{\delta }{T}_{\mathit{xxx}}^{\mathrm{res}}$, (b)$\mathit{\delta }{T}_{\mathit{xxy}}^{\mathrm{res}}$, (c) $\mathit{\delta }{T}_{\mathit{xxz}}^{\mathrm{res}}$, (d) $\mathit{\delta }{T}_{\mathit{xyz}}^{\mathrm{res}}$, (e) $\mathit{\delta }{T}_{\mathit{yyx}}^{\mathrm{res}}$, (f) $\mathit{\delta }{T}_{\mathit{yyy}}^{\mathrm{res}}$, (g) $\mathit{\delta }{T}_{\mathit{yyz}}^{\mathrm{res}}$, (h) $\mathit{\delta }{T}_{\mathit{zzx}}^{\mathrm{res}}$, (i) $\mathit{\delta }{T}_{\mathit{zzy}}^{\mathrm{res}}$, and (j) $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$.

Regarding Figure 2, the GC residuals are larger in correspondence of high frequencies in the Moho variations, e.g. subducting plates and continental oceanic boundary. A detailed future investigation needs to be performed to understand the percentage of the total signal in terms of different functionals coming from the different layers, which will reveal whether the GC functionals can observe shallower or deeper layers effects. Based on the results in Figure 2(j), the comparisons between disturbing gravitational potential and radial functionals from first-order up to third-order derivatives are discussed in Section 3.6, for providing a more regional look to better demonstrate the potential of $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ for the detection of the Moho.

3.4. Laplace parameters of the gravitational curvatures

In this section, compared to the evaluation progresses of the step-wise stripping topographic and crustal effects in Tenzer, Hamayun, and Vajda (2009a, 2009b), Tenzer, Novák, Vajda, et al. (2012), Novák and Tenzer (2013), Novák et al. (2013), and Tenzer et al. (2015), the computation precision evaluation with the Laplace parameters of the GC functionals are additionally performed during the processes of step-wise stripping crustal correction. Laplace’s equation of the GC functionals can be given as (Casotto and Fantino 2009; Šprlák and Novák 2016; Šprlák, Novák, and Pitoňák 2016; Hamáčková et al. 2016; Deng and Shen 2018a, 2018b; Deng et al. 2020):

(9) $\Delta \mathit{\delta G}{C}_1^{\chi }=\mathit{\delta }{T}_{\mathit{xxx}}^{\chi }+\mathit{\delta }{T}_{\mathit{yyx}}^{\chi }+\mathit{\delta }{T}_{\mathit{zzx}}^{\chi }$(9)

(10) $\Delta \mathit{\delta G}{C}_2^{\chi }=\mathit{\delta }{T}_{\mathit{xxy}}^{\chi }+\mathit{\delta }{T}_{\mathit{yyy}}^{\chi }+\mathit{\delta }{T}_{\mathit{zzy}}^{\chi }$(10)

(11) $\Delta \mathit{\delta G}{C}_3^{\chi }=\mathit{\delta }{T}_{\mathit{xxz}}^{\chi }+\mathit{\delta }{T}_{\mathit{yyz}}^{\chi }+\mathit{\delta }{T}_{\mathit{zzz}}^{\chi }$(11)

where $\mathit{\chi }=\mathit{res}$ means the Laplace parameters of the residual GC functionals. $\mathit{\chi }=1800$ and $\mathit{\chi }=17$ represent the Laplace parameters using the EGM2008 model with d/o 1800 and 17 in the spectral domain, respectively. $\mathit{\chi }\in${Topo/bath, Ice, USed, MSed, LSed, UCru, MCru, LCru} indicates the Laplace parameters for the different crustal effects using tesseroids in the spatial domain, where USed, MSed, and LSed are the upper, middle, and lower sediments; UCru, MCru, and LCru are the upper, middle, and lower consolidated crusts, respectively.

Theoretically, the values of the $\Delta \mathit{\delta G}{C}_1^{\chi }$, $\Delta \mathit{\delta G}{C}_2^{\chi }$, and $\Delta \mathit{\delta G}{C}_3^{\chi }$ should be zero. The precision level of the $\Delta \mathit{\delta G}{C}_1^{\mathit{res}}$, $\Delta \mathit{\delta G}{C}_2^{\mathit{res}}$, and $\Delta \mathit{\delta G}{C}_3^{\mathit{res}}$ can reveal the approximation errors of the residual GC functionals after the step-wise stripping crustal correction. The global maps of the $\Delta \mathit{\delta G}{C}_1^{\mathit{res}}$, $\Delta \mathit{\delta G}{C}_2^{\mathit{res}}$, and $\Delta \mathit{\delta G}{C}_3^{\mathit{res}}$ for the residual GC functionals are illustrated in Figure 3, and the corresponding statistical information values of the $\Delta \mathit{\delta G}{C}_1^{\mathit{res}}$, $\Delta \mathit{\delta G}{C}_2^{\mathit{res}}$, and $\Delta \mathit{\delta G}{C}_3^{\mathit{res}}$ in Figure 3 are presented in Table A4.

Figure 3. Global Laplace parameters of the (a) $\Delta \mathit{\delta G}{C}_1^{\mathit{res}}$ in Equation (9), (b) $\Delta \mathit{\delta G}{C}_2^{\mathit{res}}$ in Equation (10), and (c) $\Delta \mathit{\delta G}{C}_3^{\mathit{res}}$ in Equation (11) for the residual GC functionals.

In Table A4, the precision level of Laplace parameters of the residual GC functionals (e.g. $\Delta \mathit{\delta G}{C}_1^{\mathit{res}}$, $\Delta \mathit{\delta G}{C}_2^{\mathit{res}}$, and $\Delta \mathit{\delta G}{C}_3^{\mathit{res}}$) is about 10⁻²¹ m⁻¹s⁻². Figures 3(a–c) show that the approximation errors of the $\Delta \mathit{\delta G}{C}_1^{\mathit{res}}$, $\Delta \mathit{\delta G}{C}_2^{\mathit{res}}$, and $\Delta \mathit{\delta G}{C}_3^{\mathit{res}}$reveal the characteristics of the noise errors with the scattered points in the global map, where they are mostly about ± 0.5 x 10⁻²¹ m⁻¹s⁻².

In order to show the computation precision of the step-wise stripping crustal correction, the Laplace parameters of the GC functionals using the EGM2008 model with different d/o 1800 and 17 in the spectral domain and using tesseroid formulae for the different crustal effects in the spatial domain are shown in Table A5.

Regarding the GC Laplace parameters in the spectral domain with d/o 1800 and 17 in Table A5, the levels of the approximation errors are about 10⁻²² or 10⁻²¹ m⁻¹s⁻². For the GC Laplace parameters of the different crustal effects in the spatial domain, the approximation errors are at the levels of 10⁻²⁹, 10⁻²⁸, 10⁻²⁷ or 10⁻²⁶ m⁻¹s⁻².

Moreover, Table A4 and A5 provide the information on the quality of how Laplace’s equation is satisfied. It demonstrates that the agreement of the satisfied Laplace’s equation does not depend much on the considered mass anomalies but the used computation techniques. The precision of the Laplace parameters of the residual GC functionals (e.g. $\Delta \mathit{\delta G}{C}_1^{\mathit{res}}$, $\Delta \mathit{\delta G}{C}_2^{\mathit{res}}$, and $\Delta \mathit{\delta G}{C}_3^{\mathit{res}}$) is more influenced by the calculations with the spherical harmonic expansion of the EGM2008 model in the spectral domain (e.g. $\Delta \mathit{\delta G}{C}_1^{1800}$, $\Delta \mathit{\delta G}{C}_1^{17}$, $\Delta \mathit{\delta G}{C}_2^{1800}$, $\Delta \mathit{\delta G}{C}_2^{17}$, $\Delta \mathit{\delta G}{C}_3^{1800}$, and $\Delta \mathit{\delta G}{C}_3^{17}$) than these of the spatial domain. In order to improve the precision level of the GC Laplace parameters for the residual functionals in Table A4 from 10⁻²¹ to 10⁻²⁶ m⁻¹s⁻², the key is to improve the calculation precision with the spherical harmonic expansion in the spectral domain. In Section 3.3, the standard forward column method (Holmes and Featherstone 2002) is applied for the evaluation of the GC functionals (Hamáčková et al. 2016), where this method is accurate up to the d/o 1800. In this paper from Table A5, the precision level for the GC Laplace parameter at the d/o 1800 is about 10⁻²¹ m⁻¹s⁻², which is larger than 10⁻²² m⁻¹s⁻² of the GC Laplace parameter at the d/o 17. The truncation of d/o 1800 may be biased by aliasing. When evaluating gravity field models up to d/o 1800, the $1^{{ }^{\circ }}$ means a severe spatial under sampling. In other words, the properties of the high-degree model cannot be correctly represented by a grid of this low spatial resolution (Hirt et al. 2017). Therefore, the extended-range arithmetic approach or X-number approach up to an arbitrary degree (Fukushima 2012; Bucha and Janák 2013; Rexer and Hirt 2015) should be utilized for the evaluation of the GC functionals in the spectral domain to reach the proper computation precision.

The numerical results of the Laplace parameters of the residual GC functionals in Figure 3 and Table A4–A5 demonstrate the correctness of the progresses for the residual GC functionals during the step-wise stripping crustal correction. In a short, the used computation methodology provides a precision that is sufficient for this study.

3.5. Correlation with the Moho depth

The correlation provides general information on the suitability of each gravitational component to infer information on the Moho depth through gravity inversion. The PCC values are applied to analyze the spatial correlation between the residual gravitational functionals and Moho depths.

The PCC values ($\mathit{\gamma }$) between the residual GC functionals and the Moho depths from the global CRUST1.0 crustal model can be calculated with Equation (8)(8) $\mathit{\gamma }=\frac{{\sum }_{\mathit{i}=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{{\sum }_{\mathit{i}=1}^n{(x_i-\bar{x})}^2}\sqrt{{\sum }_{\mathit{i}=1}^n{(y_i-\bar{y})}^2}}$(8) , where the sorted numerical results from large to small are listed in Table A6. Among the PCC values of the residual GC functionals in Table A6, the $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ has the largest absolute PCC value with about 0.31.

In order to show the effects of selected d/o values for the application of spectral filtering technique on the PCC values, the d/o values (i.e. 10, 36, 45, 90, and 180) replace the d/o value 17 of the EGM2008 model in the process of the gravitational stripping correction. The absolute PCC values between the residual GC component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ and the CRUST1.0 Moho depths with different d/o values are listed in Table A7. As CRUST1.0 has not such a resolution larger than d/o 180, and by adding the d/o larger than 180, no frequency will be added. In other words, after the d/o 180 in Table A7, the correlation does not increase. The PCC value of $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ increases from 0.30 with the d/o 10 to 0.55 with the d/o 180. The change rate of the PCC value for $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ is 0.25 from the d/o 10 to 180.

In Tenzer, Hamayun, and Vajda (2009a), the absolute PCC value between the consolidated crust-stripped gravity disturbances with the CRUST2.0 Moho depths is 0.92. And the absolute PCC value between the gravity disturbances $\mathit{\delta }{g}^{\mathit{cs}}$ and the CRUST1.0 Moho depths is 0.98 in Tenzer et al. (2015). From Table A6 of this paper, the absolute PCC value between the residual GC functionals $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ and the CRUST1.0 Moho depths is 0.31 with the d/o 17, which is smaller than 0.92 (Tenzer, Hamayun, and Vajda 2009a) and 0.98 (Tenzer et al. 2015). It can be inferred that higher correlations are obtained for the first-order derivatives followed by the potential and to a much lesser extent the third-order derivatives.

3.6. Residual gravitational curvatures compared to the lower-order functionals at the Himalayan region

In order to reveal the different properties of the residual GC functionals compared to the lower-order functionals (i.e. GP, GV, and GGT), we select the Himalayan region with longitude from ${70}^{{ }^{\circ }}$ E to ${140}^{{ }^{\circ }}$E and latitude from $0^{{ }^{\circ }}$ to ${55}^{{ }^{\circ }}$N of the residual radial-radial-radial GC component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ in Section 3.3, and calculate the residual GP ($\mathit{\delta }{T}^{\mathrm{res}}$), radial GV ($\mathit{\delta }{T}_z^{\mathrm{res}}$), radial-radial GGT ($\mathit{\delta }{T}_{\mathit{zz}}^{\mathrm{res}}$) at the Himalayan region based on the same theoretical and numerical layout of this paper using the GrafLab (Bucha and Janák 2013) and Tesseroids (Uieda, Valéria, and Braitenberg 2016).

Figure 4(a) shows the Moho depths of the Himalayan region using the $1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$ CRUST1.0 model, Figure 4(b–e) depict the different gravitational functionals as $\mathit{\delta }{T}^{\mathrm{res}}$ for the GP, $\mathit{\delta }{T}_z^{\mathrm{res}}$ for the radial GV, $\mathit{\delta }{T}_{\mathit{zz}}^{\mathrm{res}}$ for the radial-radial GGT, and $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ for the radial-radial-radial GC, respectively. It should be noted that the $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ values multiply by − 1 with respect to the values in Figure 2(j) in order to keep the color consistent with other functionals. Under the same grid size of $1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$, the values of the GC component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ in Figure 4(e) can clearly reveal the Moho variations at the Himalayan region in Figure 4(a), compared to the other gravitational functionals in Figure 4(b–d). Moreover, the statistical information of the Figure 4 is presented in Table A8.

Figure 4. (a) Moho depths at the Himalayan region using the CRUST1.0 model; the calculated values of the (b) residual GP $\mathit{\delta }{T}^{\mathrm{res}}$, (c) residual GV component $\mathit{\delta }{T}_z^{\mathrm{res}}$, (d) residual GGT component $\mathit{\delta }{T}_{\mathit{zz}}^{\mathrm{res}}$, and (e) residual GC component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$.

Taking the Tibetan plateau as an example, the residual values of the GC component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ can better show the basins’ Moho depth than other residual parameters (i.e. $\mathit{\delta }{T}^{\mathrm{res}}$, $\mathit{\delta }{T}_z^{\mathrm{res}}$, and $\mathit{\delta }{T}_{\mathit{zz}}^{\mathrm{res}}$), especially for the Qaidam, Sichuan, Tarim, and Turpan basins. In other words, the radial-radial-radial GC component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ can to some extent sharpen better the detailed features of a mass anomaly, e.g. at places where the mass anomaly over the Tibetan plateau region.

4. Conclusions and outlooks

In this contribution, the crustal effects of the GC functionals and the Moho variations sensed by the GC with ETOPO1 and CRUST1.0 are investigated. The increased spatial details of the Moho variations can be detected with the $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$. The potential of the $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ can improve detailed features of mass anomalies with considerable step functions. For future studies, this could be an interesting aspect to look into it, for example, simulating the GC signal from various mass anomalies.

Particularly, the crustal effects of the disturbing GC functionals ($\mathit{\delta }{T}_{\mathit{ijk}}$) are calculated globally with ETOPO1 and CRUST1.0 models at the grid size of $1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$. Compared to the numerical methods in the spectral domain in Novák and Tenzer (2013), Novák et al. (2013) and Tenzer, Hamayun, and Vajda (2009b), Tenzer et al. (2012), and Tenzer et al. (2015), the tesseroid formulae in the spatial domain are applied in this paper because of the proper numerical efficiency and relative approximation errors at the satellite height of 260 km. The computation errors of the crustal effects of the GC functionals are evaluated with the Laplace parameters. Numerical experiments show that the precision levels for the crustal effects of the GC Laplace parameters are about 10⁻²⁹, 10⁻²⁸, 10⁻²⁷ or 10⁻²⁶ m⁻¹s⁻².

The residual gravitational functionals’ variations with the Moho depths are obtained by the gravitational stripping correction. The initial model values of the disturbing GC functionals are from the EGM2008 with truncated maximum degree 1800. After applying the spectral filtering technique, the un-modeled effects are removed by subtracting the first 17 d/o value of the model values. Meanwhile, the precision levels for the model values and un-modeled values of the EGM2008 with d/o 1800 and 17 are approximately 10 ⁻²² or 10 ⁻²¹ m ⁻¹s ⁻² for the GC Laplace parameters. Furthermore, the precision level of Laplace parameters for the residual GC functionals after the step-wise stripping crustal correction is about 10 ⁻²¹ m ⁻¹ s ⁻².

The spatial analysis with the PCC values between the residual GC functionals and the CRUST1.0 Moho depths is performed. Among the residual GC functionals, the PCC value of the residual disturbing radial component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ is largest with 0.31. According to the magnitude of the PCC values, the order of the PCC values from large to small can be divided into four levels: level 10⁻¹ for the $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$, $\mathit{\delta }{T}_{\mathit{xxz}}^{\mathrm{res}}$, and $\mathit{\delta }{T}_{\mathit{yyz}}^{\mathrm{res}}$; level 10⁻² for the $\mathit{\delta }{T}_{\mathit{yyx}}^{\mathrm{res}}$, $\mathit{\delta }{T}_{\mathit{zzx}}^{\mathrm{res}}$, and $\mathit{\delta }{T}_{\mathit{xxx}}^{\mathrm{res}}$, level 10⁻³ for the $\mathit{\delta }{T}_{\mathit{yyy}}^{\mathrm{res}}$, $\mathit{\delta }{T}_{\mathit{zzy}}^{\mathrm{res}}$, and $\mathit{\delta }{T}_{\mathit{xyz}}^{\mathrm{res}}$; and level 10⁻⁴ for the $\mathit{\delta }{T}_{\mathit{xxy}}^{\mathrm{res}}$.

Moreover, in order to investigate the un-modeled effects of the selected d/o in the spectral filtering technique, the d/o values are changed from 17 to 10, 36, 45, 90, and 180 to reveal the variations of the PCC values. With the increased d/o values, the PPC values increase as well, which means that the demanded d/o values in the spectral filtering technique for the un-modelled effects of the GC functionals should be carefully studied. There may be some shortcomings with the spectral filtering technique, which assumes that effects for d/o $\le$ 17 are solely from deeper mass anomalies, e.g. only contained in EGM2008. However, as some signals may also come from the crustal mass anomalies, rigorously, the GC effects from the crustal mass anomalies also need to be filtered in the same way as otherwise the effects are removed twice. This artifact may also explain the increased correlation when the spectral signal for the d/o > 17 is removed from EGM2008.

When calculating the crustal effects of the GC functionals, the global atmospheric effects (Anderson 1976; Tenzer, Vajda, and Hamayun 2009; Lin and Denker 2019) could be considered. The separation of topography and bathymetry, and the depth-dependent density of the seawater (Tenzer, Novák, and Gladkikh 2012; Tenzer et al. 2015) should also be included. Meanwhile, the impact of the computation height variation on the PCC values should be investigated in the future research. Moreover, the spectral characteristics of the residual GC functionals should be analyzed in ongoing research. The comparison of signals via reduction rates instead of PCC values might be more instructive. Generally, the PCC values can be high, but signal strengths could be quite different, and in this case, the PCC values are potentially misleading (Hirt et al. 2017).

The results in this study especially the data of the residual GC functionals could be applied to the gravity inversion of the global Earth’s Moho depths (Ye et al. 2016; Ye, Tenzer, and Liu 2017; Chen 2017; Uieda and Barbosa 2017; Chen, Tenzer, and Li 2018) in the next study, which could investigate whether the GC functionals help to improve the results of the Moho through the gravity inversion. Meanwhile, the higher-order gravitational potential gradients with Laplace’s equation sensed the crustal mass anomalies could also be investigated in future research (Deng et al. 2020; Deng and Ran 2021). The numerical approaches both in the spectral domain and spatial domain (e.g. X-number approach (Fukushima 2012; Bucha and Janák 2013; Rexer and Hirt 2015), Gauss-FFT method (Wu 2018a, 2018b; Bowin 2000), and conditional split and double exponential quadrature approach (Fukushima 2018; Deng 2023)) for the GC functionals could be investigated with the analysis of numerical efficiency.

Acknowledgments

We are very grateful to Prof. P. Holota, Prof. Roman Pasteka, and two anonymous reviewers for their valuable comments and suggestions, which greatly improved the manuscript. This study is supported by National Natural Science Foundation of China [Grant Nos. 42030105, 41721003, 41804012, 41631072, 41874023].

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available from the corresponding author on reasonable request.

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

 This study is supported by National Natural Science Foundation of China [Grant Nos. 42030105, 41721003, 41804012, 41631072, 41874023].

Notes on contributors

Xiao-Le Deng

Xiao-Le Deng received his PhD degree from Wuhan University in 2019. His current research mainly focuses on applications of gravity and magnetic field modeling (e.g., gravitational curvatures and magnetic curvatures) in geoscience.

Wen-Bin Shen

Wen-Bin Shen is the head and professor of Department of Geophysics, School of Geodesy and Geomatics, Wuhan University. He received his PhD degree from Graz Technical University in Austria in 1996. He holds memberships of IAG, EGU, AGU, and IUGG. His research interests focus on relativistic geodesy, Earth rotation and global change, and Earth’s free oscillation.

Michael Kuhn

Michael Kuhn is the Associate Professor of the School of Earth and Planetary Sciences, Curtin University.

Christian Hirt

Christian Hirt received his PhD degree from the University of Hannover in Germany in 2001. His research interests are modern geodesy, gravity field modeling and determination, geodetic astronomy, and atmospheric refraction.

Roland Pail

Roland Pail is the Professor of the Institute for Astronomical and Physical Geodesy, Technical University of Munich. He received his PhD degree from Graz Technical University in Austria in 1999. His research interests are physical geodesy, satellite geodesy, global and regional Earth gravity field modeling, calibration/validation of satellite observations, future gravity missions, mass transport processes in the Earth system, and synthetic Earth models.

Appendix

Table A1. Numerical values of the computation parameters in the numerical experiments.

Parameter	Notation	Magnitude	Unit
Mean Earth radius	$\mathit{R}$	6378136.3	m
Newtonian gravitational constant	$\mathit{G}$	6.67408 x 10^−11	m^3 kg^−1 s^−2
Geocentric gravitational constant	$\mathit{GM}$	3.986004415 x 10^14	m^3 s^−2
Height above reference sphere	$\mathit{h}$	260000	m
Reference crust density	${\rho }^{\mathrm{crust}}$	2670	kg m^−3
Density of the topography	${\rho }^{\mathrm{topo}}$	2670	kg m^−3
Density contrast of the bathymetry	$\mathit{\delta }{\rho }^{\mathrm{bath}}$	−1630	kg m^−3
Grid size		$1^{{ }^{\circ }}\times 1^{{ }^{\circ }}$	arc degree

Table A2. Statistical information of crustal effects of the $\mathit{\delta }{T}_{\mathit{zzz}}$ with the factor x 10⁻¹⁵ m⁻¹ s⁻² removed.

Crustal effects	Min	Max	Mean	STD
Figure 1(a) topography/bathymetry	−25.16	35.02	−0.18	5.05
Figure 1(b) ice	−11.37	5.98	−0.07	1.27
Figure 1(c) upper sediment	−3.88	2.48	−0.05	0.78
Figure 1(d) middle sediment	−3.76	2.40	−0.01	0.60
Figure 1(e) lower sediment	−1.43	0.44	0.00	0.11
Figure 1(F) upper consolidated crust	−3.08	5.42	0.01	0.76
Figure 1(g) middle consolidated crust	−4.00	4.70	0.06	0.91
Figure 1(h) lower consolidated crust	−5.10	11.59	0.21	1.32

Table A3. Statistical information of the residual GC functionals after the gravitational stripping correction with the factor × 10 ⁻¹⁵ m ⁻¹s ⁻² removed.

Parameter	Min	Max	Mean	STD
Figure 2(a) $\mathit{\delta }{T}_{\mathit{xxx}}^{\mathrm{res}}$	−120.46	88.51	0.17	8.39
Figure 2(b) $\mathit{\delta }{T}_{\mathit{xxy}}^{\mathrm{res}}$	−43.35	42.03	7.52 x 10^−11	3.52
Figure 2(c) $\mathit{\delta }{T}_{\mathit{xxz}}^{\mathrm{res}}$	−112.58	89.82	−0.01	9.23
Figure 2(d) $\mathit{\delta }{T}_{\mathit{xyz}}^{\mathrm{res}}$	−35.17	39.47	4.32 x 10^−11	3.19
Figure 2(e) $\mathit{\delta }{T}_{\mathit{yyx}}^{\mathrm{res}}$	−39.51	39.06	0.14	3.54
Figure 2(F) $\mathit{\delta }{T}_{\mathit{yyy}}^{\mathrm{res}}$	−103.89	79.39	−3.29 x 10^−10	8.26
Figure 2(g) $\mathit{\delta }{T}_{\mathit{yyz}}^{\mathrm{res}}$	−79.55	92.33	−0.14	9.10
Figure 2(h) $\mathit{\delta }{T}_{\mathit{zzx}}^{\mathrm{res}}$	−102.90	136.74	−0.31	10.64
Figure 2(i) $\mathit{\delta }{T}_{\mathit{zzy}}^{\mathrm{res}}$	−88.86	115.47	6.15 x 10^−10	10.51
Figure 2(j) $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$	−118.92	128.01	0.15	15.18

Table A4. Statistical information of the Laplace parameters for the residual GC functionals with the unit m ⁻¹s ⁻².

Parameter	Min	Max	Mean	STD
$\Delta \mathit{\delta G}{C}_1^{\mathit{res}}$	−3.90 x 10^−21	2.98 x 10^−21	−1.53 x 10^−25	1.39 x 10^−22
$\Delta \mathit{\delta G}{C}_2^{\mathit{res}}$	−1.10 x 10^−21	1.10 x 10^−21	−6.50 x 10^−25	1.39 x 10^−22
$\Delta \mathit{\delta G}{C}_3^{\mathit{res}}$	−2.03 x 10^−21	5.00 x 10^−21	1.66 x 10^−24	1.97 x 10^−22

Table A5. Statistical information of the GC Laplace parameters for the different d/o values (1800 and 17) and crustal effects with the unit m ⁻¹ s ⁻².

Parameter	Min	Max	Mean	STD
$\Delta \mathit{\delta G}{C}_1^{1800}$	−4.00 x 10^−21	3.00 x 10^−21	−5.78 x 10^−26	1.36 x 10^−22
$\Delta \mathit{\delta G}{C}_1^{17}$	−1.00 x 10^−22	1.00 x 10^−22	9.56 x 10^−26	3.16 x 10^−24
$\Delta \mathit{\delta G}{C}_1^{\mathrm{Topo}/\mathrm{bath}}$	−3.72 x 10^−27	4.40 x 10^−27	−4.87 x 10^−30	5.53 x 10^−28
$\Delta \mathit{\delta G}{C}_1^{\mathrm{Ice}}$	−1.22 x 10^−27	9.23 x 10^−28	−1.97 x 10^−30	8.21 x 10^−29
$\Delta \mathit{\delta G}{C}_1^{\mathrm{USed}}$	−3.98 x 10^−28	4.44 x 10^−28	4.06 x 10^−31	5.66 x 10^−29
$\Delta \mathit{\delta G}{C}_1^{\mathrm{MSed}}$	−3.75 x 10^−28	4.49 x 10^−28	1.37 x 10^−30	3.30 x 10^−29
$\Delta \mathit{\delta G}{C}_1^{\mathrm{LSed}}$	−1.62 x 10^−28	1.31 x 10^−28	7.75 x 10^−32	5.77 x 10^−29
$\Delta \mathit{\delta G}{C}_1^{\mathrm{UCru}}$	−4.49 x 10^−28	4.56 x 10^−28	−4.70 x 10^−31	4.32 x 10^−29
$\Delta \mathit{\delta G}{C}_1^{\mathrm{MCru}}$	−5.07 x 10^−28	5.77 x 10^−28	−1.01 x 10^−30	6.15 x 10^−29
$\Delta \mathit{\delta G}{C}_1^{\mathrm{LCru}}$	−8.80 x 10^−28	8.92 x 10^−28	4.24 x 10^−31	1.31 x 10^−28
$\Delta \mathit{\delta G}{C}_2^{1800}$	−1.00 x 10^−21	1.00 x 10^−21	−7.33 x 10^−25	1.36 x 10^−22
$\Delta \mathit{\delta G}{C}_2^{17}$	−1.00 x 10^−22	1.00 x 10^−22	−8.35 x 10^−26	3.15 x 10^−23
$\Delta \mathit{\delta G}{C}_2^{\mathrm{Topo}/\mathrm{bath}}$	−4.73 x 10^−27	6.04 x 10^−27	−1.99 x 10^−31	5.24 x 10^−28
$\Delta \mathit{\delta G}{C}_2^{\mathrm{Ice}}$	−1.09 x 10^−27	9.35 x 10^−28	−8.44 x 10^−33	7.94 x 10^−29
$\Delta \mathit{\delta G}{C}_2^{\mathrm{USed}}$	−4.92 x 10^−28	4.67 x 10^−28	−4.92 x 10^−31	5.30 x 10^−29
$\Delta \mathit{\delta G}{C}_2^{\mathrm{MSed}}$	−3.42 x 10^−28	3.50 x 10^−28	−8.37 x 10^−33	3.16 x 10^−29
$\Delta \mathit{\delta G}{C}_2^{\mathrm{LSed}}$	−9.19 x 10^−29	1.03 x 10^−28	−3.42 x 10^−33	4.90 x 10^−30
$\Delta \mathit{\delta G}{C}_2^{\mathrm{UCru}}$	−4.72 x 10^−28	4.34 x 10^−28	−3.68 x 10^−32	4.21 x 10^−29
$\Delta \mathit{\delta G}{C}_2^{\mathrm{MCru}}$	−7.00 x 10^−28	4.44 x 10^−28	−4.49 x 10^−31	5.72 x 10^−29
$\Delta \mathit{\delta G}{C}_2^{\mathrm{LCru}}$	−1.07 x 10^−27	9.03 x 10^−28	−3.17 x 10^−32	1.19 x 10^−28
$\Delta \mathit{\delta G}{C}_3^{1800}$	−2.00 x 10^−21	5.00 x 10^−21	1.46 x 10^−24	1.93 x 10^−22
$\Delta \mathit{\delta G}{C}_3^{17}$	−1.00 x 10^−22	1.00 x 10^−22	−1.94 x 10^−25	3.90 x 10^−23
$\Delta \mathit{\delta G}{C}_3^{\mathrm{Topo}/\mathrm{bath}}$	−4.23 x 10^−26	2.29 x 10^−26	5.19 x 10^−30	3.87 x 10^−27
$\Delta \mathit{\delta G}{C}_3^{\mathrm{Ice}}$	−8.05 x 10^−27	7.20 x 10^−27	4.39 x 10^−29	7.89 x 10^−28
$\Delta \mathit{\delta G}{C}_3^{\mathrm{USed}}$	−2.23 x 10^−27	3.03 x 10^−27	1.87 x 10^−29	3.82 x 10^−28
$\Delta \mathit{\delta G}{C}_3^{\mathrm{MSed}}$	−1.48 x 10^−27	2.85 x 10^−27	4.12 x 10^−30	1.82 x 10^−28
$\Delta \mathit{\delta G}{C}_3^{\mathrm{LSed}}$	−5.66 x 10^−28	5.45 x 10^−28	3.90 x 10^−31	2.22 x 10^−29
$\Delta \mathit{\delta G}{C}_3^{\mathrm{UCru}}$	−2.76 x 10^−27	1.83 x 10^−27	−1.54 x 10^−29	3.18 x 10^−28
$\Delta \mathit{\delta G}{C}_3^{\mathrm{MCru}}$	−4.11 x 10^−27	2.50 x 10^−27	−1.37 x 10^−29	5.51 x 10^−28
$\Delta \mathit{\delta G}{C}_3^{\mathrm{LCru}}$	−7.33 x 10^−27	4.35 x 10^−27	−2.98 x 10^−29	1.22 x 10^−28

Table A6. The sorted absolute PCC values ($\mathit{\gamma }$) from large to small between the Moho depths from CRUST1.0 model and the residual GC functionals after the gravitational stripping correction.

Parameter	$\mathit{\gamma }$
$\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$	0.31
$\mathit{\delta }{T}_{\mathit{xxz}}^{\mathrm{res}}$	0.26
$\mathit{\delta }{T}_{\mathit{yyz}}^{\mathrm{res}}$	0.25
$\mathit{\delta }{T}_{\mathit{yyx}}^{\mathrm{res}}$	1.80 x 10^−2
$\mathit{\delta }{T}_{\mathit{zzx}}^{\mathrm{res}}$	1.48 x 10^−2
$\mathit{\delta }{T}_{\mathit{xxx}}^{\mathrm{res}}$	1.11 x 10^−2
$\mathit{\delta }{T}_{\mathit{yyy}}^{\mathrm{res}}$	7.63 x 10^−3
$\mathit{\delta }{T}_{\mathit{zzy}}^{\mathrm{res}}$	6.16 x 10^−3
$\mathit{\delta }{T}_{\mathit{xyz}}^{\mathrm{res}}$	3.18 x 10^−3
$\mathit{\delta }{T}_{\mathit{xxy}}^{\mathrm{res}}$	4.60 x 10^−4

Table A7. The absolute PCC values ($\mathit{\gamma }$) between the residual GC component $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$ and CRUST1.0 Moho depth with the different d/o for the spectral filtering technique in the process of the gravitational stripping correction.

d/o	$\mathit{\gamma }$
10	0.30
17	0.31
36	0.31
45	0.33
90	0.47
180	0.55

Table A8. Statistical information of the residual functionals at the Himalayan region with the factor x 10⁻¹⁵ m⁻¹ s⁻² removed for the $\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$.

Parameter	Min	Max	Mean	STD
Moho depths	9.05 km	74.81 km	32.45 km	14.20 km
$\mathit{\delta }{T}^{\mathrm{res}}$	−7017.62 m^2 s^−2	6789.05 m^2 s^−2	750.19 m^2 s^−2	3903.05 m^2 s^−2
$\mathit{\delta }{T}_z^{\mathrm{res}}$	−263.69 mGal	654.07 mGal	147.06 mGal	221.22 mGal
$\mathit{\delta }{T}_{\mathit{zz}}^{\mathrm{res}}$	−4.12 Eötvös	9.90 Eötvös	1.16 Eötvös	2.60 Eötvös
$\mathit{\delta }{T}_{\mathit{zzz}}^{\mathrm{res}}$	−125.42	118.92	1.45	23.14