On Moho Determination by the Vening Meinesz-Moritz Technique

Abstract

This chapter describes a theory and application of satellite gravity and altimetry data for determining Moho constituents (i.e. Moho depth and density contrast) with support from a seismic Moho model in a least-squares adjustment. It presents and applies the Vening Meinesz-Moritz gravimetric-isostatic model in recovering the global Moho features. Internal and external uncertainty estimates are also determined. Special emphasis is devoted to presenting methods for eliminating the so-called non-isostatic effects, i.e. the gravimetric signals from the Earth both below the crust and from partly unknown density variations in the crust and effects due to delayed Glacial Isostatic Adjustment as well as for capturing Moho features not related with isostatic balance. The global means of the computed Moho depths and density contrasts are 23.8±0.05 km and 340.5 ± 0.37 kg/m3, respectively. The two Moho features vary between 7.6 and 70.3 km as well as between 21.0 and 650.0 kg/m3. Validation checks were performed for our modeled crustal depths using a recently published seismic model, yielding an RMS difference of 4 km.

1. Introduction

Traditionally, the structure of the Earth’s interior is divided according to its chemical and physical properties into crust, mantle, outer core and inner core. The oceanic crust ranges from 5 to 10 km depth, while the continental crust ranges from 35 to 70 km depth. The layer below the crust is the mantle, which is the thickest layer of the Earth. It can be divided into the upper (extending down to 660 km from the Earth’s surface) and lower mantle (down to 2900 km beneath the surface). The innermost layer of the Earth is the core, which can be decomposed into the outer and inner core. (A modern decomposition of the Earth’s interior is based on its main mechanical properties: the lithosphere and asthenosphere, of main interest in global geodynamics, plate tectonics and motion, but not for this study).

The geoscientist typically uses three sources of information to figure out the interior of the Earth’s structure:

The first source is understood by direct evidence from rock samples by drilling projects. In this way, the scientist attempts to drill holes in the Earth’s surface, to a maximum depth of about 12 km, and explode rocks for inferring the conditions within the Earth’s interior. The drilling method is severely limited, because it is difficult to drill a deep hole due to the high pressure and temperature, and it is also a very time-consuming and expensive technology (see [1]).

The second source includes the records of seismic waves, which are generated, for example, by earthquakes, explosions, volcanoes and other natural sources. Accordingly, specialists can detect information about the Earth’s interior, e.g. depth to density discontinuities, through detailed analysis of seismic data. Also, by studying the velocity of the wave, it can to some extent be used for estimating the density of the medium. At this point it deserves to be mentioned that the seismic data are also expensive to collect and therefore sparse and in-homogeneously distributed around the Earth (see [2]).

The third set of information in modeling the Earth’s interior is the recent gravity field models, generated through modern satellite gravity missions such as Challenging Mini-satellite Payload (CHAMP), Gravity Recovery and Climate Experiment (GRACE) and Gravity field and steady state Ocean Circulation Explorer (GOCE), which can provide a global and homogeneous coverage of data. An improvement can also be obtained in the accuracy and spatial resolution of these models by combining them with airborne and ground-based gravity data as well as satellite altimetry data over the oceans. Other important sources for studying Earth’s interior are its magnetic field and meteorites.

1.1 Background of Moho modeling

The primary interface of the Earth’s interior is the boundary between the Earth’s crust and mantle, which is called the Mohorovičić discontinuity (or Moho). This discontinuity was first discovered in 1909 by the Croatian seismologist Andrija Mohorovičić, when analyzing seismograph records of an earthquake in the Kapula valley, namely P-waves (compressional waves) and S-waves (shear waves). He noticed that the P-waves, which travel deeper into the Earth, moved faster than those that travel nearer the surface. Accordingly, he concluded that the Earth is not homogeneous, and at a specific depth there must be a boundary surface, which distinguishes two media with different compositions, and by which the seismic waves propagate with different velocities (see [3]).

Currently the Moho interface can be studied using two main methods: the gravimetric and seismic ones. These methods cannot provide exactly the same results, as they are based on different hypotheses, different types, qualities and spatial distributions of data (see, e.g. [4, 5]).

The seismic methods are the major traditional techniques in modeling the thickness of the Earth’s crust (the Moho depth, MD), where the base of the crust is defined as the Moho. Another Moho constituent is the Moho Density Contrast (MDC), which can be estimated from the change of velocity of a seismic wave passing through the Moho boundary. Models based on seismic data can be locally very accurate but useless in areas without adequate seismic observations, particularly over large portions of the oceans. In addition, the seismic data acquisition is costly with lack of global coverage [6].

In contrast, while using satellite gravity data, information on the Moho can be inferred from a uniform and global data set. However, Moho models based on gravity data are in general characterized by simplified hypotheses to guarantee the uniqueness of the solution of the inverse gravitational problem (see, e.g. [7]). As we will show in Section 2.1, gravity data alone cannot separate the MD from the MDC, but additional information is needed to solve this problem. In any case, due to the complementary information described above, a combined gravimetric-seismic method could be fruitful in modeling the Moho.

Much research using seismic surveys for recovering the Moho interface has been performed in the last decades. For instance, [8, 9] compiled global Moho models based on seismic data analysis, and [10] estimated the MD using seismic surface waves. For global studies the most frequently used crustal models are the CRUST2.0 [11] and CRUST1.0 models [12], compiled with 2° × 2° and 1° × 1° resolutions, respectively. More recently, [13] developed a global crustal thickness model and velocity structure from geostatistical analysis of seismic data, and we hereafter call this model CRUST19.

Over large areas of the world with a sparse coverage of seismic data, in particular at sea, a gravimetric-isostatic or combined gravimetric/seismic method can be prosperous. For example, [14] modified the Airy/Heiskanen theory ([15], Section 3.4) by introducing a regional isostatic compensation model based on a thin plate lithospheric flexure model [16, p. 114]. [17, Section 8] generalized the Vening Meinesz hypothesis from a regional to global compensation. [7] expressed the Vening Meinesz-Moritz (VMM) problem as that of solving a non-linear Fredholm integral equation, and presented some solutions for recovering the MD. The VMM method was also followed up by some additional theoretical studies, such as methods for estimating the MDC [6] and for reducing the Bouguer gravity anomaly for non-isostatic effects [18, 19]. [20] demonstrated that the MD estimated from the isostatic gravity disturbance based on solving the VMM model has a better agreement with the CRUST2.0 seismic model than those computed by the isostatic gravity anomaly. Their argument was also theoretically explained by [21]. [22] estimated the MD and MDC using a combination of the CRUST2.0 and a GOCE global gravity models. [23] showed that the application of the Bouguer gravity disturbance and the no-topography correction in the VMM model to determine the MD provides very similar results, suggesting the preference of the gravity disturbance to the traditional Bouguer gravity anomaly for gravity inversion. [4, 5] computed combined Moho constituent model according to the VMM method. [24] estimated a new MDC model named MDC2018, using the marine gravity field from satellite altimetry in combination with a seismic-based crustal model and Earth’s topographic/bathymetric data. Finally, [25] estimated a combined Moho model for marine areas via satellite altimetric - gravity and seismic crustal models.

1.2 Gravimetric-isostatic Moho models

Isostasy is an important concept in Earth sciences describing the state of equilibrium (or mass balance) to which the mantle tends to balance the mass of the crust in the absence of external disturbing forces. “When a certain area of the crust reaches the state of isostasy, it is said to be in isostatic equilibrium (or balance), and the depth at which isostatic equilibrium prevails is called the depth of compensation” [26].

However, the transport of material over the Earth’s surface, such as glaciers, volcanism, and sedimentation, etc., are factors that disturb isostasy, yielding so-called non-isostatic effects (NIEs).

Four principle models of isostasy related with the crustal depth and/or density can briefly be listed as those of (a) Airy/Heiskanen (A/H; [27, 28, 29]), (b) Pratt/Hayford (P/H; [30, 31]), (c) Vening Meinesz (VM; [14]), and (d) the Vening Meinesz-Moritz (VMM; [7, 17]). Common for the isostatic models is that the Bouguer gravity anomaly $\mathrm{\Delta }{g}^b$ (or disturbance $\delta g^b$) is fully compensated by a compensation attraction below the crust such that the isostatic gravity anomaly and disturbance vanish:

$$\mathrm{\Delta }{g}^I=\delta g^I=0.$$E1

A/H and P/H are local models, implying that the compensation attraction operates along the vertical of the observation point, implying that the sum of the masses of the crust and its compensation along each vertical is assumed to be constant from place to place.

The A/H model assumes a constant crustal density, and variations in topographic height is compensated by variations in the depth of the crust. That is, the mass excess of topography is compensated by the mass deficit of mountain roots in the upper mantle. In ocean areas anti-roots of mantle material compensates for the light mass of the ocean.

The P/H model assumes a constant depth of compensation of the solid Earth topography (including negative topography over oceans), while the density of the topography varies with topographic height.

Due to the elasticity of the Earth’s crust these local models are not very realistic. Hence, [14] modified the A/H model by introducing a model with a regional compensation in which mass loads and unloads are balanced by a gentle bending or flexure of the crust over a regional area. [17] generalized the VM model from a regional to a global compensation with a spherical sea level approximation. [7] and, finally, [6] generalized the VMM model to allow for variations both in crustal density and depth. In this way the VMM can be seen as a generalization of both the A/H and P/H models with global isostatic compensations by variations of both mountain root and crustal density.

Below we will present the least-squares theory for determining a combined VMM-seismic model for both MD and MDC. The theory is finally applied in a new global model.

2. The VMM theory

2.1 Solution for the product of Moho depth and density contrast $\left(D\mathrm{\Delta }\rho \right)$

In H. Moritz’ original publication [17] the problem is to determine the MD $\left(D\right)$ such that the compensation attraction ($A_C$) fully compensates the Bouguer gravity anomaly. Here we employ this condition in the last part of Eq. (1), which can be written (cf. [7, 21])

$$\delta g^I=\delta g^B+A_C=0,$$E2

where $\delta g^B$ is the Bouguer gravity disturbance (i.e., the free-air gravity disturbance after removal of the topographic attraction).

The VMM technique uses both gravimetric and seismic data in a least squares combination to determine the MD (D) and/or MDC $\left(\mathrm{\Delta }\rho \right)$. The method assumes that the crust is in isostatic balance, implying that the isostatic gravity anomaly $\left(\mathrm{\Delta }{g}^I\right)$ and disturbance $\left(\delta g^I\right)$ vanish at each point on the Earth’s surface as in Eq. (1) above. Note that the compensation attraction is a function of both MD and MDC. Approximating the Earth’s surface by a sphere of radius R, one obtains after several manipulations of Eq. (2) the following equation in D for a constant Δρ:

$$\mathit{RG}\mathrm{\Delta }\rho \underset{\sigma }{\iint }K\left(\psi ,s\right)\mathit{d\sigma }=f,$$E3

where G is the gravitational constant, $K\left(\psi ,s\right)$ is an integral kernel function with arguments $\psi =$ geocentric angle between integration and computation points and $s=1-D/R$, and $f=-\left(\delta g^b+A_{C0}\right)/G$. Here $A_{C0}$ is zero-degree harmonic of the compensation attraction (which does not affect the Moho undulation). Eq. (3) is a non-linear Fredholm integral equation of the first kind, which has the following first- and second-order solutions:

$$D_1=\frac{1}{4\mathit{\pi ∆p}}\sum _{n=0}^{\infty }\frac{2n+1}{n+1}\sum _{m=-n}^n{f}_{\mathit{nm}}{Y}_{\mathit{nm}}$$E4

where $Y_{\mathit{nm}}$ is a fully-normalized spherical harmonic, $f_{\mathit{nm}}$ is the corresponding coefficient given by the Bouguer gravity disturbance f, and

$${\left(D_2\right)}_P={\left(D_1\right)}_P+\frac{{\left(D_1^2\right)}_P}{R}-\frac{1}{32\mathit{R\pi }}\underset{\sigma }{\iint }\left[\frac{{\left(D_1^2\right)}_Q-{\left(D_1^2\right)}_P}{{sin}^3{\psi }_{\mathit{PQ}}}\right]d{\sigma }_Q.$$E5

Here subscripts P and Q denote computation and integration points, respectively, $f_{\mathit{nm}}$ is the spherical harmonic coefficient of f. Note that the integral contributes significantly only locally around the computation point. The formula can be improved by a few steps of iteration:

$$D_P^{k+1}=D_P^k+\frac{{\left(D_P^k\right)}^2}{R}-\frac{1}{32\mathit{R\pi }}\underset{\sigma }{\iint }\left[\frac{{\left(D_Q^k\right)}^2-{\left(D_P^k\right)}^2}{{sin}^3{\psi }_{\mathit{PQ}}}\right]d{\sigma }_Q;k=0,1,2,\dots ,$$E6

where $D_P^0=D_1$ at point P determined by Eq. (4).

As the isostatic balance of the crust is hardly valid for crustal blocks of diameter smaller than, say, 100 km ([32], p.195), the upper limit of the series in Eq. (4) of should not exceed $n_2$ =180. Also, as we shall see later, the low-degree harmonics in D₁, say, below $n_1$ =10, are not contributing to the isostatic balance but are due to mass anomalies in the Earth’s interior below the crust.

The integrals in Eqs. (5) and (6) are local, as the integrand quickly vanishes with distance away from the computation point. Hence a flat earth approximation may be relevant (See [6]).

If the MDC varies laterally, the following 2nd-order approximation of Eq. (3) can be found in the spectral domain (cf. [6]) when introducing the notation $\chi =D\mathrm{\Delta }\rho$

$$f_{\mathit{nm}}=4\pi \frac{n+1}{2n+1}\left[{\chi }_{\mathit{nm}}+\frac{n+2}{2R}{\left(\mathit{\chi D}\right)}_{\mathit{nm}}\right]$$E7

and, after summing up, one obtains:

$$\chi =\sum _{n=n_1}^{n_2}\left[\frac{2n+1}{n+1}\frac{f_n}{4\pi }-\frac{n+2}{2R}{\left(\mathit{\chi D}\right)}_n\right],$$E8

where $f_n$ and ${\left(\mathit{\chi D}\right)}_n$ are the Laplace harmonics

$$\left(\begin{array}{c}{f}_n\\ {\left(\mathit{\chi D}\right)}_{\mathit{n}}\end{array}\right)=\sum _{m=-n}^n\left(\begin{array}{c}{f}_{\mathit{nm}}\\ {\left(\mathit{\chi D}\right)}_{\mathit{nm}}\end{array}\right)Y_{\mathit{nm}}.$$E9

Using the approximation

$$\mathit{\chi D}\approx {\chi }^2/\mathrm{\Delta }{\rho }_P,$$E10

one obtains from (8) the iterative formula

$${\chi }_P^{k+1}=\sum _{n=n_1}^{n_2}\left[\frac{2n+1}{n+1}\frac{f_n}{4\pi }-\sum _{m=-n}^n\frac{n+2}{2R\mathrm{\Delta }{\rho }_P}{\left({\chi }^2\right)}_n^k\right],;k=0,1,2,\dots$$E11

where ${\chi }_P^0$ is the first-order solution:

$${\chi }^0=\frac{1}{4\pi }\sum _{n=n_1}^{n_2}\frac{2n+1}{n+1}{f}_n.$$E12

Alternatively, we may present Eq. (11) by the iterative formula:

$${\chi }_P^{k+1}={\chi }_P^0+I_P^k,;k=0,1,2,\dots ,$$E13

where

$$I_P^k=\frac{1}{R\mathrm{\Delta }{\rho }_P}\underset{\sigma }{\iint }\frac{{\left({\chi }_P^k\right)}^2-{\left({\chi }_Q^k\right)}^2}{{sin}^3\psi }d{\sigma }_Q.$$E14

Again, this integral is very local, which suggests the use of a flat-Earth approximation. Also, assuming that $\left(n_2+2\right)D_0/\left(2R\right)<1$, Eq. (7) leads to the approximate solution:

$${\chi }_P\approx \frac{1}{4\pi }\sum _{n=n_1}^{n_2}\frac{2n+1}{n+1}\frac{f_n}{1+\left(n+2\right)D_0/\left(2R\right)}.$$E15

Note that the solution ${\chi }_P$ is the product of the MD and MDC. If one of the parameters is known, the other can be determined by the equation. Hence, gravity data alone cannot be used to distinguish between the two Moho constituents. Hence, additional information, e.g., from seismic and/or geological data, is needed to separate the two. However, as we shall see later, usually such data is not taken for granted in the VMM technique, but the gravity data used in Eq. (8) is typically applied to improve a priori Moho constituents in a least-squares procedure.

The solution (8) can be derived from Eq. (1), and from the inversion of a 3-D Newton integral. See Appendix A.

2.2 A least-squares solution for both the Moho depth and the Moho density contrast

The Moho component $\chi$, the product of D and Δρ, can be estimated from Eq. (11) or (12) and applied as an observation together with seismic data for solving both the MD and the MDC in a least-squares adjustment. Then a linear set of equations including gravimetric data ($l_1$), and seismic data for MD ($l_2$) and MDC ($l_3$) can be written for each pixel (P):

$$\mathrm{\Delta }{\rho }_P\mathit{dD}+D_{P}d\mathrm{\Delta }\rho =l_1-{\epsilon }_1$$E16

$$\mathit{dD}=l_2-{\epsilon }_2$$E17

$$d\mathrm{\Delta }\rho =l_3-{\epsilon }_3,$$E18

Where dD and $d\mathrm{\Delta }\rho$ are the (unknown) corrections to the initial values $D_P$ and $\mathrm{\Delta }{\rho }_P$, $l_1=\chi -{\chi }_P-I_P^k$, where $I_P^0=0$, and ${\epsilon }_i$ are the errors of the observations. In matrix form the adjustment system can be written

$$\mathbf{A}\mathbf{X}=\mathbf{L}-\mathbf{\epsilon },$$E19

where

$$\mathbf{A}=\left[\begin{array}{cc}\mathrm{\Delta }{\rho }_P& D_P\\ 1& 0\\ 0& 1\end{array}\right],\mathbf{X}=\left[\begin{array}{c}\mathit{dD}\\ d\mathrm{\Delta }\rho \end{array}\right]\text{and}\mathbf{L}=\left[\begin{array}{c}{l}_1\\ l_2\\ l_3\end{array}\right].$$E20

Assuming that the observation errors are random with expectation zero and covariance matrix Q, the weighted least squares solution of this system becomes:

$$\widehat{\mathbf{X}}\mathbf{=}{\left({\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{A}\right)}^{\mathbf{-}\mathbf{1}}{\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}^{\mathbf{-}\mathbf{1}}\mathbf{L}.$$E21

From this result, the adjusted MD and MDC for point P are obtained by:

$$\widehat{D}=D_P+d\widehat{D}\text{and}\mathrm{\Delta }\widehat{\rho }=\mathrm{\Delta }{\rho }_P+d\mathrm{\Delta }\widehat{\rho }.$$E22

As the first equation $\left(l_1\right)$ is a linearization, it could make sense to iterate the adjustment procedure by replacing the previous initial values $D_P$ and $\mathrm{\Delta }{\rho }_P$ in Eq. (16) by their adjusted values $\widehat{D}$, $\mathrm{\Delta }\widehat{\rho }$ and repeat the above computation procedure until sufficient convergence.

3. Uncertainty estimations

First, the result of the least-squares procedure depends on the quality and weighting of the gravity and seismic observations. The weights should be selected as proportional to the inverse standard errors (STEs) of the observations squared. The STEs of seismic data is, hopefully, provided along with the data files. For the gravity data we derive the global mean STE in Section 3.1. In Section 3.2 we propagate the data errors to error estimates in the VMM least-squares results of Moho constituents. Finally, in Section 3.3 a method for validating the modeled Moho undulations is presented.

3.1 The uncertainty in the gravimetric-isostatic observation equation

Assuming that there are no systematic errors and disregarding 2nd –order terms in Eq. (8), one obtains the error in $\chi$ by simple error propagation from Eq. (12):

$${\epsilon }_{\chi }=\frac{1}{4\pi }\sum _{n=n_1}^{n_2}\frac{2n+1}{n+1}{\mathit{df}}_n,$$E23

where ${\mathit{df}}_n$ is the error in $f_n$. Then it follows that the global Root Mean Square Error (RMSE) of $\chi$ becomes

$$\mathit{\text{RMSE}}\left(\chi \right)=\frac{1}{4\pi }E\left\{\underset{\sigma }{\iint }\left({\epsilon }_{\chi }^2\right)\mathit{d\sigma }\right\}=\frac{1}{4\pi }\sqrt{\sum _{n_1}^{n_2}{\left(\frac{2n+1}{n+1}\right)}^2{\mathit{dc}}_n},$$E24

where $E\left\{\right\}$ denotes the statistical expectation of the term in the bracket, and ${\mathit{dc}}_n$ are the error degree variances of the gravity disturbances. Using this formula with harmonics between 10 and 180 of the XGM2019e gravity field model (see [33]), the RMSE value becomes 1.17 × 10⁴ kg/m².

3.2 The uncertainties in VMM Moho depth and density contrast

Assuming that all observation errors are stochastic with expectation zero, an error propagation of the least squares solution in Eq. (21) yields that the covariance matrix of $\widehat{\mathbf{X}}$ becomes

$${\mathbf{Q}}_{\mathbf{XX}}={\sigma }_0^2{\left({\mathbf{A}}^{\mathrm{T}}{\mathbf{Q}}^{-1}\mathbf{A}\right)}^{-1}$$E25

where ${\sigma }_0^2$ is the variance of unit weight, which can be unbiasedly estimated by

$$s_0^2={\mathbf{L}}^{\mathrm{T}}{\mathbf{Q}}^{-1}\left(\mathbf{L}\mathbf{-}\mathbf{{\rm A}}\widehat{\mathbf{X}}\right).$$E26

Note that there is no denominator in Eq. (26), because in the present adjustment example with 3 observations and 2 unknowns per pixel there is only 1 degree of freedom.

3.3 Verification of the solutions

First, we will find an estimate of the variance ${\sigma }_x^2$ of the solution x for the MD or MDC by assuming that we know another solution y with variance ${\sigma }_y^2$ . If both solutions have vanishing expected errors, the solution becomes

$${\sigma }_x^2={\sigma }_y^2+E\left\{x^2-y^2\right\}.$$E27

The correlation coefficient between x and y follows from

$$k=\left[{\sigma }_x^2+{\sigma }_y^2-E\left\{{\left(x-y\right)}^2\right\}\right]/\left(2{\sigma }_x{\sigma }_y\right)$$E28

One can also plot the t-test parameter of the normalized (and unitless) difference between x and y:

$$T=\frac{x-y}{\sqrt{{\sigma }_x^2+{\sigma }_y^2-2k{\sigma }_x{\sigma }_y}}$$E29

to study the expected difference.

To verify Eqs. (27)–(29), one may start from the substitutions that the true value for x and y is given by

$$\overline{x}=x-e_x=y-e_y,$$E30

where $e_x$ and $e_y$ are random errors with zero-expectations.

In practice, x and y are the Moho quantities at a pixel estimated from two models, and the expectation operator should be replaced by the (weighted) mean value over the central and surrounding pixels. Note that the solution in Eq. (27) is independent on whether x and y are correlated or not. Eq. (29) can be used in an approximate t-test to judge whether the estimates x and y from the two models are statistically equal or not, if they are (weighted) mean values.

4. Corrections to gravimetric data

Nowadays, the Earth’s gravity field has been recognized as an important source of information about the Earth’s structure. Such data contain both short- and long-wavelength features, i.e., signals from the topographic and bathymetry geometries and density heterogeneities in the topography, ice caps, sediment basins and also in the mantle and core/mantle topography variations.

The long-wavelength contribution to the gravity field, say to spherical harmonic degree and order 10, may be assumed to be related to the mantle and below located heterogeneities.

To isolate the gravity data caused only by the geometry and density contrast of the Moho interface, all aforementioned signal contributors to the gravity data must be removed by applying the so-called stripping corrections and NIEs [34] and NIEs (see section 4.2). Another gravity correction corresponds to the gravimetric effect of filling-up all oceans with masses to a standard density of 2670 kg/m³. Finally, by removing also normal gravity from the resulting stripped free-air gravity observation, one obtains the refined Bouguer gravity disturbance. As a result, the ideal stripped Bouguer gravity disturbance can be explained as caused by a spherical Earth without solid Earth topography and mass anomalies below the crust.

4.1 Crustal density corrections

In order to compute the stripped refined Bouguer gravity disturbance, i.e. free-air gravity disturbance corrected for topography, bathymetry, ice thickness and sediment basins (i.e. stripping corrections), [34] developed and applied a uniform mathematical formalism of computing the gravity corrections of the density variations within the Earth’s crust. This operation can be summarized as the correction

$$\delta g^{\mathit{\text{TBIS}}}=\delta g^t+\delta g^b+\delta g^i+\delta g^s$$E31

where $\delta g^t$ is the topographic gravity correction, and $\delta g^b$, $\delta g^i$ and $\delta g^s$ are the stripping gravity corrections due to the ocean (bathymetry), ice and sediment density variations, respectively.

Applying a spherical approximation of the Earth, each gravity correction on the right-hand side of Eq. (31) can be computed using the following spherical harmonic series:

$$\delta g^q\left(P\right)=\frac{\mathit{GM}}{R^2}\sum _{n=0}^{n_{\max }}\left(n+1\right)\sum _{m=-n}^n{c}_{\mathit{nm}}^q{Y}_{\mathit{nm}}\left(\mathrm{P}\right),$$E32

with superscript q being one of t, b, i or s, and $\mathit{GM}=3986005\times {10}^8{\mathrm{m}}^3{\mathrm{s}}^{-2}$ is the geocentric gravitational constant. The coefficient $c_{\mathit{nm}}^q$ of a particular volumetric mass density (or density contrast) layer q (i.e., topography, bathymetry, glacial ice and sediments) is defined by:

$$c_{\mathit{nm}}^q=\frac{2}{\left(2n+1\right)}\frac{1}{{\rho }_e}\left[\frac{{\left({\rho }^q{L}_q\right)}_{\mathit{nm}}}{R}+\frac{n+2}{2}\frac{{\left({\rho }^q{L}_q^2\right)}_{\mathit{nm}}}{R^2}+\dots \right]$$E33

where ${\rho }^q$ is the Earth’s mean mass density, and the coefficients (${\rho }^q{L}_q^i$) are evaluated (from discrete data of density ${\rho }^q$ and thickness $L_q$) by applying a discretization to the following integral convolution

$${\left({\rho }^q{L}_q\right)}_{\mathit{nm}}=\frac{1}{4\pi }\underset{\sigma }{\iint }{\rho }^q{L}_q^j{Y}_{\mathrm{n}m}\mathit{d\sigma },j=1,2,\dots ,n$$E34

4.2 Non-isostatic effects

It is important to remind the reader that in general the crust is not in complete isostatic equilibrium, and the observed gravity data are not only generated by the topographic/isostatic masses, but also from those in the deep Earth interior, that leads to non-isostatic effects (NIEs) (see [18, 19, 35]).

According to [7], the major part of the long-wavelengths of the geopotential undulation is caused by density variations in the Earth’s mantle and core/mantle topography variations. Such NIEs could be the contribution of different factors, such as crustal thickening/thinning, thermal expansion of mass of the mantle [36], Glacial Isostatic Adjustment (GIA), plate flexure ([16], p. 114), and effect of other phenomena. This implies that this contribution to gravity will lead to systematic errors/NIEs of the computed Moho topography. Hence the NIEs should also be corrected on the isostatic gravity disturbance.

Assuming that the seismic Moho model CRUST1.0 is known and correct, the gravity effect of the NIEs can be determined by:

$$\delta g^{\mathit{NIE}}=\frac{\mathit{GM}}{R^2}\sum _{n=0}^{n_{\max }}\left(n+1\right)\sum _{m=-n}^n{c}_{\mathit{nm}}^{\mathit{NIE}}{Y}_{\mathit{nm}}\left(\mathrm{P}\right)$$E35

where

$$c_{\mathit{nm}}^{\mathit{NIE}}=c_{\mathit{nm}}^{\text{CRUST}1.0}-c_{\mathit{nm}}^{\mathit{VMM}}$$E36

Here $c_{\mathit{nm}}^{\mathit{NIE}},c_{\mathit{nm}}^{\mathit{VMM}},c_{\mathit{nm}}^{\text{CRUST}1.0}$ are the spherical harmonic coefficients of the gravity disturbances of the NIE, VMM and CRUST1.0, respectively.

The isostatic equilibrium equation in Eq. (2) is then rewritten as:

$$\delta g^I\left(P\right)=\delta g_B^{\mathit{\text{TBISN}}}\left(P\right)+A_C\left(P\right)=0.$$E37

Here $\delta g_B^{\mathit{\text{TBISN}}}$ is the refined Bouguer gravity disturbance corrected for the gravitational contributions of topography and density variations of the oceans, ice, sediments and NIEs, i.e. by Eq. (31).

4.3 Glacial isostatic adjustment (GIA)

Delayed GIA (DGIA) expresses the delayed adjustment process of the Earth to an equilibrium state when former ice sheet loads have vanished. The ongoing adjustment of the Earth’s body to the redistribution of ice and water masses is evident in various phenomena, which have been studied to infer the extent and amount of the former ice masses, to reconstruct the sea level during a glacial cycle and to constrain rheological properties of the Earth’s interior. Here we aim at answering the question whether the effect of the gravimetric DGIA correction is significant for Moho determination in Fennoscandia. Usually, this effect is part of and reduced by the general NIE correction, but one may also estimate the DGIA effect on gravity as a separate correction by the harmonic window:

$$\delta g^{\mathit{\text{DGIA}}}=\gamma \sum _{n=10}^{23}\left(n+1\right)\sum _{m=-n}^n{A}_{\mathit{nm}}{Y}_{\mathit{nm}},$$E38

where $\gamma$ is normal gravity, $Y_{\mathit{nm}}$ and $A_{\mathit{nm}}$ are spherical harmonics and coefficients of the gravitational potential (see [37]). Here the limits of the series are based on the optimum correlation between the present land uplift and the gravity field in the region.

5. A global VMM solution

The main gravimetric input data to be used in the following VMM Moho model is the global Earth gravitational field model (e.g. XGM2019e) in the harmonic window from $n_1=10$ to $n_2=180$. The gravity disturbance data were corrected for the gravitational signals of mass density variations due in different layers of the Earth’s crust (i.e. stripping gravity corrections) and for the gravity contribution from deeper masses below the crust (i.e. non-isostatic effects). The NIEs were computed using the seismic crustal model CRUST1.0, and the stripping corrections for different crustal heterogeneous data utilized the global topographic models DTM2006 and Earth2014. The preliminary gravimetric Moho solution was combined with the CRUST1.0 model in a least-squares procedure (see Section 2.2). The adjustment was performed globally for each $1^{\circ }\times 1^{\circ }$-block.

Quantities $\mathit{\delta g}$ (mGal)	Max	Mean	Min	STD
$\mathit{\delta g}$	285.85	−0.44	−281.40	23.84
$\delta g^t$	255.13	−71.06	−647.61	105.98
$\delta g^b$	721.60	332.91	110.28	165.02
$\delta g^I$	325.78	21.84	−2.61	56.57
$\delta g^S$	185.31	45.48	−0.02	32.47
$\delta g^{\mathit{NIE}}$	248.70	−134.65	−497.98	69.98
$\delta g^{\mathit{\text{TBISN}}}$	562.82	128.87	−620.54	178.10

Table 1.

Statistics of global estimates of the gravity disturbances, stripping gravity corrections and NIEs. STD is the standard deviation of the estimated quantity over the blocks. $\mathit{\delta g}$ is the gravity disturbance computed by the XGM2019e coefficients. $\delta g^t$, $\delta g^b$, $\delta g^I$ and $\delta g^S$ are the topographic/bathymetric, ice and sediment stripping gravity corrections derived from the CRUST1.0, respectively. $\delta g^{\mathit{NIE}}$ is the non-isostatic effect. $\delta g^{\mathit{\text{TBISN}}}$ is the refined Bouguer gravity disturbance after applying the topographic and stripping gravity corrections due to the ocean, ice and sediment density variations.

The statistics of the stripping gravity corrections and refined Bouguer gravity disturbance are presented in Table 1. It shows the largest corrections for bathymetry and NIE, but also ice cap corrections have some extreme values. The sum of the corrections varies roughly within $\pm 600$ mGal with the STD of 178 mGal.

Figure 1 depicts the Bouguer gravity disturbances corrected for the ocean (bathymetry), ice, sediment variations and the NIEs, respectively. As one can see from the figure, these features can drastically change the Bouguer gravity disturbance from the free-air disturbance over oceans due to the application of the bathymetric stripping gravity correction. It also changes in central Greenland and Antarctica due to the applied ice density variation stripping gravity correction (Figure 1d). In Figure 1e one can see large stripping corrections in sediment basins, and the NIEs are also very significant (Figure 1f).

Figure 1.

(a) The free-air gravity disturbance computed using the XGM2019e coefficients complete to degree 180 of spherical harmonics, (b) the topographic gravity correction, (c) the bathymetric stripping gravity correction, (d) the ice density variation stripping gravity correction, (e) the sediments density variation stripping gravity corrections, (f) non-isostatic effects and (g) refined Bouguer gravity disturbances after applying the above corrections. (h) the DGIA effect in Fennoscandia. Unit: mGal.

Figure 1g shows the refined Bouguer gravity disturbance after applying the above corrections. This disturbance has a span of about $\pm$ 500 mGal, to be compared with the approximate span of $\pm$ 250 mGal of the free-air disturbance. Notable is the large positive disturbances on the oceans corresponding to the effect of filling the oceans with topographic masses. The DGIA effect, demonstrated for Fennoscandia and, depicted in Figure 1h, is very small compared to other corrections.

In the least-squares procedure of the combined VMM solution the weights of the two types of data were chosen as follows. The weights of the gravity disturbances were estimated from their inverse variances by Eq. (23), while the weights for CRUST1.0 data were those published in [12]. Figures 2 and 3 depict the results of the MD and density contrast undulations and their estimated standard errors. Their extreme values for continental and oceanic crusts and mean values are reported in Table 2.

Figure 2.

(a) The MD estimated from combined approach, and (b) its standard error. (unit km).

Figure 3.

(a) The MDC estimated by combined approach, and (b) its standard error. (unit kg/m³).

Quantities		Max.	Mean	Min.	STD
MD (km)	Global	70.26	23.78	7.55	13.17
	Ocean	43.19	14.98	7.55
	Land	70.26	40.03	18.37
STE MD (km)	Global	8.15	1.20	0.05	0.94
	Ocean	7.34	2.06	0.05
	Land	8.15	2.49	1.05
MDC (kg/m^3)	Global	649.99	340.49	20.98	100.90
	Ocean	637.36	281.01	20.98
	Land	649.99	440.01	69.34
STE MDC (kg/m^3)	Global	132.26	17.44	0.09	14.17
	Ocean	99.98	35.21	0.09
	Land	132.26	38.65	19.09

Table 2.

Statistics of global estimates of MD and MDC in the VMM approach for 1° × 1° block data. STD is the standard deviation. STE is the standard error obtained in the least-squares adjustment. Units for MD and MDC are km and kg/m³, respectively.

To validate the STE of the VMM solution for crustal depth, we determined the global mean of it by Eq. (27) using the seismic model CRUST19. The result is 1.73 km, which is in fair agreement with the 1.20 km given in Table 2. Also, as one can see from Figure 4 the test parameter in Eq. (29) for validating the VMM solution of MD from the seismic model CRUST19 is mainly in the range $\pm$ 1, which suggest rather close agreements of estimated MDs and their error estimates. (Note that $E\left\{{\mathrm{T}}^2\right\}=1$, implies that assumed variance components are correct and the expected MDs of the two models are the same).

Figure 4.

Validation of the VMM MD solution by Eq. (29) and CRUST19 model. (the scale is unitless).

6. Discussion and final remarks

The study of the Moho discontinuity has been a crucial topic in inferring the dynamics of the Earth’s interior for a long time. In general, the Moho can be studied with profitable results through seismic data. However, due to the sparsity of seismic data in parts of the world, it has not been well determined. With the advent of satellite missions, it has been possible to recover the Moho constituents via satellite gravity observations based on an isostatic model.

So far, various isostatic models have been presented for recovering the Moho constituents, but it was not clarified which one is most appropriate to employ for geophysical and geodynamical purposes. The preliminary and simplest isostatic models proposed are the classical ones with local or regional compensation. However, those models cannot realistically image the actual Moho undulation. This is because they assume a uniform crustal density, disregarding the density irregularities distributed within the crust and sub-crust. Understanding this important role of Moho recovery has been in the center of the discussions by many geoscientists during the last decades.

Here we have determined the Moho constituents and their uncertainties based on the VMM technique using both gravimetric and seismic data on a global scale to a resolution of 1° × 1°. The combination of the gravimetric and seismic data in one approach as well as the joint adjustment of MD and density contrast are expected to significantly improve the total result.

The basic VMM method is based on the hypothesis that the isostatic gravity disturbance vanish. However, this is the case only if the gravity component is reduced such that there are no signals from the Earth’s interior below the crust. The major problem in this reduction is therefore to distinguish and remove those signals, which we utilize by estimating and removing the NIEs with the help from CRUST1.0 seismic model.

The second step is to combine the gravimetric data, propagated in the VMM technique to a linear equation (with MD and MDC as the unknowns), with a seismic model, CRUST1.0. This is performed by a weighted least-squares adjustment, block by block, which has the advantage that the standard error of the unknowns can also be estimated block-wise. The weights of the gravity disturbances were based on the error estimates by Eq. (23), while the weights for CRUST1.0 data were those published in [12].

Our estimated results can be summarized as follows. The global means of MD and MDC are 23.8 $\pm$ 0.05 km and 340.5 $\pm$ 0.37 kg/m³, respectively, ranging between 7.6–70.3 km and 21.0–650.0 kg/m³. The MD results were validated by the recent CRUST19 seismic model, showing that the differences between the models vary within the extremes −23.4 and 32.9 km, with a global average of 0.91 km and an RMS fit of 4 km. The normalized differences were generally within the limits $\pm 1$, which should be regarded as acceptable.

Acknowledgments

This study was supported by project no. 187/18 of the Swedish National Space Agency (SNSA).

Appendix A

Let us assume that the compensation attraction in Eq. (2) is generated by a density contrast $\mathrm{\Delta }\rho$ between the constant reference depth $D_0$ and the actual depth D. Assuming that the density contrast may change only laterally, it follows from the Newton integral in 3D, that the compensation potential becomes:

$${\left(V,c\right)}_P=G\underset{\sigma }{\iint }\mathrm{\Delta }\rho \underset{R-D}{\overset{R-D_0}{\int }}\frac{r^2\mathit{dr}}{l_P}\mathit{d\sigma }=G\underset{\sigma }{\iint }\mathrm{\Delta }\rho \underset{R-D}{\overset{R}{\int }}\frac{r^2\mathit{dr}}{l_P}\mathit{d\sigma }+G\underset{\sigma }{\iint }\mathrm{\Delta }\rho \underset{R}{\overset{R-D_0}{\int }}\frac{r^2\mathit{dr}}{l_P}\mathit{d\sigma },$$EA.1

where the last integral term is a constant, global mean value. Disregarding this term (which does not contribute to the Moho undulation) the integral can be written in the spectral domain after integration with respect to r and setting $r_P=R$ (sea level radius):

$${\left(V,c\right)}_P={\mathit{GR}}^2\sum _{n=0}^{\infty }\frac{1}{n+3}\underset{\sigma }{\iint }\frac{\mathrm{\Delta }\rho }{n+3}\left[1-{\left(1-\frac{D}{R}\right)}^{n+3}\right]P_n\left(cos\psi \right)\mathit{d\sigma }$$EA.2

Considering the addition theorem of fully normalized spherical harmonics (Heiskanen and Moritz 1967, p. 33):

$$P_n\left(cos{\psi }_{\mathit{PQ}}\right)=\frac{1}{2n+1}\sum _{m=-n}^n{Y}_{\mathit{nm}}\left(P\right)Y_{\mathit{nm}}\left(Q\right),$$

one obtains

$$\delta T_P={\mathit{GR}}^2\sum _{n=0}^{\infty }\frac{1}{\left(2n+1\right)\left(n+3\right)}\sum _{m=-n}^n{Y}_{\mathit{nm}}\left(P\right)\underset{\sigma }{\iint }\mathrm{\Delta }\rho \left[1-{\left(1-\frac{D}{R}\right)}^{n+3}\right]Y_{\mathit{nm}}\mathit{d\sigma }.$$EA.3

As D is small vs. R, one may expand the last bracket in this equation a la Taylor as

$$\frac{1}{n+3}\left[1-{\left(1-\frac{D}{R}\right)}^{n+3}\right]=\frac{D}{R}-\frac{n+2}{2}{\left(\frac{D}{R}\right)}^2+\frac{\left(n+2\right)\left(n+1\right)}{6}{\left(\frac{D}{R}\right)}^3+\dots ,$$EA.4

and by inserting this series in Eq. (A.3) one obtains after integration

$$\begin{array}{l}{\left(V,c\right)}_P=4\mathit{\pi GR}\sum _{n=0}^{\infty }\frac{1}{2n+1}\times \\ \sum _{m=-n}^n\left[{\left(D\mathrm{\Delta }\rho \right)}_{\mathit{nm}}-\frac{n+2}{2R}{\left(D^2\mathrm{\Delta }\rho \right)}_{\mathit{nm}}+\frac{\left(n+2\right)\left(n+1\right)}{6R^2}{\left(D^3\mathrm{\Delta }\rho \right)}_{\mathit{nm}}-\dots \right]Y_{\mathit{nm}}\left(P\right)\end{array}$$EA.5

where ${\left(\right)}_{\mathit{nm}}$ are spherical harmonic coefficients. As the compensation potential coefficients are related to those of the compensation attraction Ac by

$${\left(V,c\right)}_{\mathit{nm}}=\frac{R}{n+1}{\left(\mathit{Ac}\right)}_{\mathit{nm}},$$EA.6

one obtains the spectral equation from (A.5):

$$\frac{R}{n+1}{\left(\mathit{Ac}\right)}_{\mathit{nm}}=\frac{4\mathit{\pi GR}}{2n+1}\left[{\left(D\mathrm{\Delta }\rho \right)}_{\mathit{nm}}-\frac{n+2}{2R}{\left(D^2\mathrm{\Delta }\rho \right)}_{\mathit{nm}}+\frac{\left(n+2\right)\left(n+1\right)}{6R^2}{\left(D^3\mathrm{\Delta }\rho \right)}_{\mathit{nm}}-\dots \right]$$EA.7

By comparing the spectra of both sides and summing up all harmonics and considering Eq. (2),

$$\mathit{D}\mathrm{\Delta }\rho =\frac{1}{4\mathit{\pi }}\sum _{n=0}^{\infty }\frac{2n+1}{n+1}\sum _{m=-n}^n{f}_{\mathit{nm}}{Y}_{\mathit{nm}}+\mathrm{\Delta },$$EA.8

where f_nm = −δ $g_{\mathit{nm}}^B$/G and

$$\mathrm{\Delta }=\sum _{n=0}^{\infty }\sum _{m=-\mathrm{n}}^{\infty }\left[\frac{n+2}{2R}{\left(D^2\mathrm{\Delta }\rho \right)}_{\mathit{nm}}-\frac{\left(n+2\right)\left(n+1\right)}{6R^2}{\left(D^3\mathrm{\Delta }\rho \right)}_{\mathit{nm}}+\mathrm{..}\right]Y_{\mathit{nm}}$$EA.9

accounts for higher order terms in the series.

In practical application for Moho feature determination the lower limit for the summation is found to be about 10 (as the lower harmonics are related with deep Earth gravity anomalies) and the upper limit should not exceed about 180.