Open access peer-reviewed chapter

On Moho Determination by the Vening Meinesz-Moritz Technique

Written By

Lars Erik Sjöberg and Majid Abrehdary

Submitted: October 26th, 2020 Reviewed: March 26th, 2021 Published: May 7th, 2021

DOI: 10.5772/intechopen.97449

Chapter metrics overview

308 Chapter Downloads

View Full Metrics


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.


  • crustal depth
  • Moho density contrast
  • Moho depth
  • Vening Meinesz-Moritz method

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 Δgb (or disturbance δgb) is fully compensated by a compensation attraction below the crust such that the isostatic gravity anomaly and disturbance vanish:


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 DΔρ

In H. Moritz’ original publication [17] the problem is to determine the MD D such that the compensation attraction (AC) 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])


where δgB 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 Δρ. The method assumes that the crust is in isostatic balance, implying that the isostatic gravity anomaly ΔgI and disturbance δgI 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 Δρ:


where G is the gravitational constant, Kψs is an integral kernel function with arguments ψ= geocentric angle between integration and computation points and s=1D/R, and f=δgb+AC0/G. Here AC0 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:


where Ynm is a fully-normalized spherical harmonic, fnm is the corresponding coefficient given by the Bouguer gravity disturbance f, and


Here subscripts P and Q denote computation and integration points, respectively, fnm 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:


where DP0=D1 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 n2 =180. Also, as we shall see later, the low-degree harmonics in D1, say, below n1 =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 χ=DΔρ


and, after summing up, one obtains:


where fn and χDn are the Laplace harmonics


Using the approximation


one obtains from (8) the iterative formula


where χP0 is the first-order solution:


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




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


Note that the solution χ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 χ, 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 (l1), and seismic data for MD (l2) and MDC (l3) can be written for each pixel (P):


Where dD and dΔρ are the (unknown) corrections to the initial values DP and ΔρP, l1=χχPIPk, where IP0=0, and εi are the errors of the observations. In matrix form the adjustment system can be written




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


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


As the first equation l1 is a linearization, it could make sense to iterate the adjustment procedure by replacing the previous initial values DP and ΔρP in Eq. (16) by their adjusted values D̂, Δρ̂ 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 χ by simple error propagation from Eq. (12):


where dfn is the error in fn. Then it follows that the global Root Mean Square Error (RMSE) of χ becomes


where E denotes the statistical expectation of the term in the bracket, and dcn 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 × 104 kg/m2.

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 X̂ becomes


where σ02 is the variance of unit weight, which can be unbiasedly estimated by


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 σx2 of the solution x for the MD or MDC by assuming that we know another solution y with variance σy2 . If both solutions have vanishing expected errors, the solution becomes


The correlation coefficient between x and y follows from


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


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


where ex and ey 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/m3. 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


where δgt is the topographic gravity correction, and δgb, δgi and δgs 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:


with superscript q being one of t, b, i or s, and GM=3986005×108m3s2 is the geocentric gravitational constant. The coefficient cnmq of a particular volumetric mass density (or density contrast) layer q (i.e., topography, bathymetry, glacial ice and sediments) is defined by:


where ρq is the Earth’s mean mass density, and the coefficients (ρqLqi) are evaluated (from discrete data of density ρq and thickness Lq) by applying a discretization to the following integral convolution


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:




Here cnmNIE,cnmVMM,cnmCRUST1.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:


Here δgBTBISN 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:


where γ is normal gravity, Ynm and Anm 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 n1=10 to n2=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×1-block.

Quantities δg (mGal)MaxMeanMinSTD

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. δg is the gravity disturbance computed by the XGM2019e coefficients. δgt, δgb, δgI and δgS are the topographic/bathymetric, ice and sediment stripping gravity corrections derived from the CRUST1.0, respectively. δgNIE is the non-isostatic effect. δgTBISN 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 ±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 ± 500 mGal, to be compared with the approximate span of ± 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/m3).

MD (km)Global70.2623.787.5513.17
STE MD (km)Global8.
MDC (kg/m3)Global649.99340.4920.98100.90
STE MDC (kg/m3)Global132.2617.440.0914.17

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/m3, 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 ± 1, which suggest rather close agreements of estimated MDs and their error estimates. (Note that ET2=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 ± 0.05 km and 340.5 ± 0.37 kg/m3, respectively, ranging between 7.6–70.3 km and 21.0–650.0 kg/m3. 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 ±1, which should be regarded as acceptable.



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 Δρ between the constant reference depth D0 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:


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 rP=R (sea level radius):


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


one obtains


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


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


where nm are spherical harmonic coefficients. As the compensation potential coefficients are related to those of the compensation attraction Ac by


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


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


where fnm = −δ gnmB/G and


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.


  1. 1. Kearey P, Brooks M, Hill I. An introduction to geophysical exploration. John Wiley & Sons; 2002 Apr 26
  2. 2. Anderson DL. Theory of the Earth. Blackwell scientific publications; 1989 Jan 1
  3. 3. Abrehdary M. Recovering Moho parameters using gravimetric and seismic data. KTH (Doctoral dissertation, Ph. D. thesis, TRITA SoM 2016–02, 56 pp., Stockholm, Sweden)
  4. 4. Abrehdary M, Sjöberg LE, Bagherbandi M. Combined Moho parameters determination using CRUST1. 0 and Vening Meinesz-Moritz model. Journal of Earth Science. 2015 Aug;26(4):607–16
  5. 5. Abrehdary M, Sjöberg LE, Bagherbandi M, Sampietro D. Towards the Moho depth and Moho density contrast along with their uncertainties from seismic and satellite gravity observations. Journal of Applied Geodesy. 2017 Dec 1;11(4):231–47
  6. 6. Sjöberg LE, Bagherbandi M. A method of estimating the Moho density contrast with a tentative application of EGM08 and CRUST2. 0. Acta Geophysica. 2011 Jun;59(3):502–25
  7. 7. Sjöberg LE. Solving Vening Meinesz-Moritz inverse problem in isostasy. Geophysical Journal International. 2009 Dec 1;179(3):1527–36
  8. 8. Meier U, Curtis A, Trampert J. Global crustal thickness from neural network inversion of surface wave data. Geophysical Journal International. 2007 May 1;169(2):706–22
  9. 9. Shapiro NM, Ritzwoller MH. Monte-Carlo inversion for a global shear-velocity model of the crust and upper mantle. Geophysical Journal International. 2002 Oct 1;151(1):88–105
  10. 10. Lebedev S, Adam JM, Meier T. Mapping the Moho with seismic surface waves: a review, resolution analysis, and recommended inversion strategies. Tectonophysics. 2013 Dec 8;609:377–94
  11. 11. Bassin C. The current limits of resolution for surface wave tomography in North America. EOS Trans. AGU. 81: Fall Meet. Suppl., Abstract. 2000
  12. 12. Laske G, Masters G, Ma Z, Pasyanos M. Update on CRUST1. 0—A 1-degree global model of Earth’s crust. InGeophys. res. abstr 2013 Apr (Vol. 15, p. 2658)
  13. 13. Szwillus W, Afonso JC, Ebbing J, Mooney WD. Global crustal thickness and velocity structure from geostatistical analysis of seismic data. Journal of Geophysical Research: Solid Earth. 2019 Feb;124(2):1626–52
  14. 14. Meinesz FV. Une nouvelle methode pour la reduction isostatique regionale de l’intensite de la pesanteur. Bulletin geodesique. 1931 Jan 1;29(1):33–51
  15. 15. Heiskanen WA, Moritz H, Physical Geodesy WH. Freeman and Company. San Francisco. 1967
  16. 16. Watts AB. Isostasy and Flexure of the Lithosphere Cambridge Univ. Press, Cambridge. 2001
  17. 17. Moritz H. The figure of the Earth: theoretical geodesy and the Earth’s interior. Karlsruhe: Wichmann. 1990
  18. 18. Bagherbandi M, Sjöberg LE. Non-isostatic effects on crustal thickness: A study using CRUST2. 0 in Fennoscandia. Physics of the Earth and Planetary Interiors. 2012 Jun 1;200:37–44
  19. 19. Bagherbandi M, Sjöberg LE. Improving gravimetric–isostatic models of crustal depth by correcting for non-isostatic effects and using CRUST2. 0. Earth-Science Reviews. 2013 Feb 1;117:29–39
  20. 20. Tenzer R, Bagherbandi M. Reformulation of the Vening-Meinesz Moritz inverse problem of isostasy for isostatic gravity disturbances. International Journal of Geosciences. 2012;3(5A):918–29
  21. 21. Sjöberg LE. On the isostatic gravity anomaly and disturbance and their applications to Vening Meinesz–Moritz gravimetric inverse problem. Geophysical Journal International. 2013 Jun 1;193(3):1277–82
  22. 22. Reguzzoni M, Sampietro D, Sansò F. Global Moho from the combination of the CRUST2. 0 model and GOCE data. Geophysical Journal International. 2013 Oct 1;195(1):222–37
  23. 23. Sjöberg LE, Bagherbandi M, Tenzer R. On gravity inversion by no-topography and rigorous isostatic gravity anomalies. Pure and Applied Geophysics. 2015 Oct;172(10):2669–80
  24. 24. Abrehdary M, Sjöberg LE, Sampietro D. Contribution of satellite altimetry in modelling Moho density contrast in oceanic areas. Journal of Applied Geodesy. 2019 Jan 28;13(1):33–40
  25. 25. Abrehdary M, Sjöberg LE. Estimating a combined Moho model for marine areas via satellite altimetric-gravity and seismic crustal models. Studia Geophysica et Geodaetica. 2020 Jan;64(1):1–25
  26. 26. Sjöberg, L. E. and Bagherbandi, M. Isostasy–Geodesy. In Encyclopedia of Geodesy. Springer International Publishing Switzerland 2014
  27. 27. Airy GB. III. On the computation of the effect of the attraction of mountain-masses, as disturbing the apparent astronomical latitude of stations in geodetic surveys. Philosophical Transactions of the Royal Society of London. 1855 Dec 31(145):101–4
  28. 28. Heiskanen WA. Untersuchungen über Schwerkraft und Isostasie. Buchdruckerei-aktiengesellschaft Sana; 1924
  29. 29. Heiskanen WA, Airy SG. New Isostaic Tables for the Reduction of Gravity Values Calculated on the Basis of Airy’s Hypothesis. Suomalainen Tideakatemia; 1938
  30. 30. Hayford JF. Geodesy: The Figure of the Earth and Isostasy from Measurements in the United States. US Government Printing Office; 1909
  31. 31. Pratt JH. I. On the attraction of the Himalaya Mountains, and of the elevated regions beyond them, upon the plumb-line in India. Philosophical Transactions of the Royal Society of London. 1855 Dec 31(145):53–100
  32. 32. Turcotte DL, Schubert G. Geodynamics. Cambridge university press; 2002 Mar 25
  33. 33. Zingerle P, Pail R, Gruber T, Oikonomidou X. The combined global gravity field model XGM2019e. Journal of Geodesy. 2020 Jul;94(7):1-2
  34. 34. Tenzer R, Chen W, Tsoulis D, Bagherbandi M, Sjöberg LE, Novák P, Jin S. Analysis of the refined CRUST1. 0 crustal model and its gravity field. Surveys in geophysics. 2015 Jan;36(1):139–65
  35. 35. Abrehdary M, Sjöberg LE, Bagherbandi M. The spherical terrain correction and its effect on the gravimetric-isostatic Moho determination. Geophysical Journal International. 2016 Jan 1;204(1):262–73
  36. 36. Kaban MK, Schwintzer P, Reigber C. A new isostatic model of the lithosphere and gravity field. Journal of Geodesy. 2004 Dec 1;78(6):368–85
  37. 37. Abrehdary M, Sjöberg LE. A New Moho Depth Model for Fennoscandia with Special Correction for the Glacial Isostatic Effect. Pure and Applied Geophysics. 2021 Feb 15:1-2

Written By

Lars Erik Sjöberg and Majid Abrehdary

Submitted: October 26th, 2020 Reviewed: March 26th, 2021 Published: May 7th, 2021