Moho Density Contrast in Central Eurasia from GOCE Gravity Gradients

Abstract

Seismic data are primarily used in studies of the Earth’s inner structure. Since large parts of the world are not yet sufficiently covered by seismic surveys, products from the Earth’s satellite observation systems have more often been used for this purpose in recent years. In this study we use the gravity-gradient data derived from the Gravity field and steady-state Ocean Circulation Explorer (GOCE), the elevation data from the Shuttle Radar Topography Mission (SRTM) and other global datasets to determine the Moho density contrast at the study area which comprises most of the Eurasian plate (including parts of surrounding continental and oceanic tectonic plates). A regional Moho recovery is realized by solving the Vening Meinesz-Moritz’s (VMM) inverse problem of isostasy and a seismic crustal model is applied to constrain the gravimetric solution. Our results reveal that the Moho density contrast reaches minima along the mid-oceanic rift zones and maxima under the continental crust. This spatial pattern closely agrees with that seen in the CRUST1.0 seismic crustal model as well as in the KTH1.0 gravimetric-seismic Moho model. However, these results differ considerably from some previously published gravimetric studies. In particular, we demonstrate that there is no significant spatial correlation between the Moho density contrast and Moho deepening under major orogens of Himalaya and Tibet. In fact, the Moho density contrast under most of the continental crustal structure is typically much more uniform.

1. Introduction

The Earth’s satellite observation systems have been used in various geoscience studies of the Earth’s interior and processes. The global elevation models derived from processing the Shuttle Radar Topography Mission (SRTM) data are used, for instance, to compute the topographic gravity correction in context of the gravimetric interpretation of sedimentary basins as well as other subsurface density structures and/or density interfaces. Apart from climatic studies, the satellite-altimetry observations are also used to determine the marine gravity data. Since the marine gravity field is highly spatially correlated with the ocean-floor relief at a certain wavelength-band, these gravity data are used to predict the ocean-floor depths [1,2,3]. By analogy with applying the topographic gravity correction, the ocean-floor depths are used to compute the bathymetric-stripping gravity correction. The satellite-gravity observations have also been used to interpret the Earth’s inner structure. This becomes particularly possible with the advent of three dedicated satellite-gravity missions, namely the Challenging Mini-satellite Payload (CHAMP) [4,5,6], the Gravity Recovery and Climate Experiment (GRACE) [7] and the Gravity field and steady-state Ocean Circulation Explorer (GOCE) [8,9]. The latest gravitational models derived from these satellite missions have a spatial resolution about 66 to 80 km (in terms of a half-wavelength). Moreover, these gravitational models have (almost) global and homogeneous coverage with well-defined stochastic properties.

Methods for a gravimetric Moho recovery from topographic, bathymetric, gravitational and crustal structure models have been developed and applied by a number of authors (e.g., [10,11,12,13,14,15,16,17]). In these studies, a uniform Moho density contrast has often been assumed. The results of seismic studies, however, revealed that the Moho density contrast varies significantly. The continental Moho density contrast 200 kg·m⁻³ in Canada was reported by [18]. Regional seismic studies [19,20] based on using the wave-receiver functions indicate that the density contrast regionally varies as much as from 160 kg·m⁻³ (for the mafic lower crust) to 440 kg·m⁻³ (for the felsic lower crust), with an apparently typical value 440 kg·m⁻³ for the craton. The Moho density contrast information is also included in the global crustal model CRUST1.0 [21]. According to this seismic model, the Moho density contrast (computed as the density difference between the upper mantle and the lower crust) globally varies between 10 and 610 kg·m⁻³, while this range is between 340 and 790 kg·m⁻³ when computed with respect to the reference crustal density 2670 kg·m⁻³ (as often used in gravimetric studies).

An attempt to incorporate the variable Moho density contrast in gravimetric methods for a determination of the Moho depth has been done by some authors. [22] took into consideration the lateral as well as radial density variation within the crystalline crustal layers and simultaneously adjusted all the densities and estimated the Moho depth. A method for a simultaneous determination of the Moho depth and density contrast developed by [23] was based on solving the Vening Meinesz-Moritz’s (VMM) inverse problem of isostasy [24,25,26,27]. This method was applied to estimate the Moho parameters globally and regionally (e.g., [28]). It was also demonstrated that the gravimetric determination of the Moho depth is more accurate when using the variable Moho density contrast [29]. The gravimetric results confirmed large variations of the Moho density contrast. According to [23], the Moho density contrast varies globally from 81.5 kg·m⁻³ (in the Pacific region) to 988 kg·m⁻³ (beneath the Tibetan Plateau). A similar range of values between 82 and 965 kg·m⁻³ was reported by [30]. It was also shown [31] that the Moho density contrast under the oceanic crust is highly spatially correlated with the ocean-floor age.

The gravimetrically-determined Moho density contrast [23,30] differs significantly from the CRUST1.0 values especially under major orogens in central Asia. To examine this aspect, we investigated the Moho density contrast at the study area which comprises most of the Eurasian tectonic plate and includes also surrounding oceanic and continental plates. For this purpose, we developed and applied a novel approach which utilizes a least-squares technique for solving the condition equations based on minimizing residuals between the gravimetric and seismic Moho parameters, particularly specified for a product of the Moho depth and density contrast.

2. Method

A functional relation between the refined gravity data and the Moho density contrast is defined here by means of solving the VMM inverse problem of isostasy. We note that this functional relation was already derived by [26], but using a slightly different approach than that presented here. We further extend this definition for finding the Moho density contrast from the in-orbit GOCE gravity-gradient data. We then propose a least-squares technique for solving the VMM problem based on combining the vertical gravity gradients and seismic model and finally investigate a spatial behavior of the integral kernel used for a regional gravity-gradient data inversion.

2.1. Moho Density Contrast from Gravity Disturbance

The VMM isostatic problem was formulated in the following generic form [26,27,32]

$$-GR{\displaystyle \underset{\sigma }{\iint }\Delta \rho K\left(\psi ,T\right)d\sigma =\delta g^i}$$

(1)

where G is the Newton’s gravitational constant, R is the Earth’s mean radius, $\Delta \rho$ is the (variable) Moho density contrast, T is the Moho depth, $\delta g^i$ is the isostatic gravity disturbance, σ is the unit sphere, and $d\sigma$ is the surface integration element. The integral kernel K in Equation (1) is a function of the spherical distance ψ and the Moho depth T. Its spectral representation reads [26]

$$K\left(\psi ,T\right)={\displaystyle \sum _{n=0}^{\infty }\frac{n+1}{n+3}\left[1-{\left(1-\frac{T}{R}\right)}^{n+3}\right]P_n\left(\text{cos} \psi \right)}$$

(2)

where P_n is the Legendre polynomial of degree n.

The isostatic gravity disturbances $\delta g^i$ on the right-hand side of Equation (1) are computed from the gravity disturbances $\delta g$ by subtracting the gravitational contributions of topography $g^{\mathrm{T}}$, bathymetry $g^{\mathrm{B}}$ and sediments $g^{\mathrm{S}}$ and adding the compensation attraction $A_{C0}$. Hence

$$\delta g^i=\delta g-g^{\mathrm{T}}-g^{\mathrm{B}}-g^{\mathrm{S}}+A_{C0}$$

(3)

The compensation attraction $A_{C0}$ in Equation (3) is computed from [26]

$$A_{C0}=\frac{4\pi G\Delta {\rho }_{0}R}{3}\left[{\left(1-\frac{T_0}{R}\right)}^3-1\right]\approx -4\pi G\Delta {\rho }_0{T}_0$$

(4)

where $T_0$ and $\Delta {\rho }_0$ are mean values of the Moho depth and density contrast respectively.

To solve the VMM problem for finding the Moho density contrast $\Delta \rho$, we first simplified the integral kernel K in Equation (2). By applying a Taylor series (up to a first-order term) to ${\left(1-T/R\right)}^{n+3}$, we get

$$K\left(\psi ,T\right)=\frac{T}{R}{\displaystyle \sum _{n=0}^{\infty }\left(n+1\right)P_n\left(\text{cos} \psi \right)}$$

(5)

Substitution from Equation (5) to Equation (1) then yields

$$-G{\displaystyle \sum _{n=0}^{\infty }\left(n+1\right)}{\displaystyle \underset{\sigma }{\iint }\left(\Delta \rho T\right)P_n\left(\text{cos} \psi \right)d\sigma =\delta g^i}$$

(6)

The surface integral on the left-hand side of Equation (6) is defined in terms of the Laplace harmonic ${\left(\Delta \rho T\right)}_n$ as follows [33]

$${\displaystyle \underset{\sigma }{\iint }\left(\Delta \rho T\right)P_n\left(\text{cos} \psi \right)d\sigma =}\frac{4\pi }{2n+1}{\left(\Delta \rho T\right)}_n$$

(7)

Inserting from Equation (7) to Equation (6), we obtained the following relation

$$-4\pi G{\displaystyle \sum _{n=0}^{\infty }\frac{n+1}{2n+1}}{\left(\Delta \rho T\right)}_n=\delta g^i$$

(8)

The coefficients ${\left(\Delta \rho T\right)}_n$ in Equation (8) are generated from the Laplace harmonics $\delta g_n^i$ of the isostatic gravity field as follows

$${\left(\Delta \rho T\right)}_n=-\frac{A_{C0}}{4\pi G}{\delta }_{n0}-\frac{1}{4\pi G}\frac{2n+1}{n+1}\left(\delta g_n-g_n^{\mathrm{T}}-g_n^{\mathrm{B}}-g_n^{\mathrm{S}}\right)$$

(9)

where $g_n^{\mathrm{T}}$, $g_n^{\mathrm{B}}$ and $g_n^{\mathrm{S}}$ are, respectively, the Laplace harmonics of the topographic, bathymetric and sediment gravitational contributions, and ${\delta }_{n0}$ is Kronecker’s delta.

The expression in Equation (9) is finally rearranged into the following form

$$\Delta \rho T=\Delta {\rho }_0{T}_0-\frac{1}{4\pi G}{\displaystyle \sum _{n=0}^{\infty }\frac{2n+1}{n+1}}\left(\delta g_n-g_n^{\mathrm{T}}-g_n^{\mathrm{B}}-g_n^{\mathrm{S}}\right)$$

(10)

As seen in Equation (10), the gravity disturbances and gravity corrections are computed in the spectral domain, while the compensation attraction $A_{C0}$ is computed approximately according to Equation (4) from the a priori values of the Moho depth and density contrast.

2.2. Moho Density Contrast from Gravity Gradient

The functional model for finding the Moho density contrast in Equation (10) is now redefined for the vertical gravity gradient (i.e., the second-order radial derivative of the disturbing potential $V_{rr}$). It is worth mentioning that results (not presented here) revealed that the contribution of horizontal gravity-gradient components on the Moho parameters is completely negligible.

The relation between the Laplace harmonics $V_{rr,n}$ and $\delta g_n$ reads

$$\delta g_n=\frac{r^2}{R}\frac{1}{\left(n+2\right)}{\left(\frac{r}{R}\right)}^{n+1}{V}_{rr,n}$$

(11)

where r is the geocentric radius. The application of a scaling factor r/R in Equation (11) allows solving the VMM problem for the gravity-gradient data at an arbitrary point on or above the geoid. Note that the generic expression in Equation (1) defines the VMM problem only on the geoid surface (approximated by a sphere of radius R). The downward continuation of gravity data observed at the topographic surface (or satellite altitudes) onto the geoid is then required in prior of solving the VMM problem.

Substitution from Equation (11) to Equation (10) yields

$$\Delta \rho T=\Delta {\rho }_0{T}_0-\frac{1}{4\pi G}\left[{\displaystyle \sum _{n=0}^{\infty }\frac{2n+1}{n+1}}\left(-g_n^{\mathrm{T}}-g_n^{\mathrm{B}}-g_n^{\mathrm{S}}\right)+W\right]$$

(12)

where the parameter W is defined by

$$W=\frac{r^2}{R}{\displaystyle \sum _{n=0}^{\infty }\frac{2n+1}{\left(n+1\right)\left(n+2\right)}}{\left(\frac{r}{R}\right)}^{n+1}{V}_{rr,n}$$

(13)

The Laplace harmonics $W_n$ of the parameter W are generated from the respective harmonics $V_{rr,n}$ of the vertical gravity gradient $V_{rr}$ using the following expression

$$W_{n}R\frac{\left(n+1\right)\left(n+2\right)}{\left(2n+1\right)}{\left(\frac{R}{r}\right)}^{n+1}=r^2{V}_{rr,n}$$

(14)

Alternatively, the Laplace harmonics $W_n$ can be defined in the following form [33]

$$W_n=\frac{2n+1}{4\pi }{\displaystyle \underset{\sigma }{\iint }W{P}_n\left(\text{cos} \psi \right)d\sigma }$$

(15)

Substituting from Equation (15) back to Equation (14) and applying further simplifications, we arrived at

$$\frac{R}{4\pi }{\displaystyle \underset{\sigma }{\iint }WL\left(r,\psi \right)d\sigma }=r^2{V}_{rr}$$

(16)

where the integral kernel L reads

$$L\left(r,\psi \right)={\displaystyle \sum _{n=0}^{\infty }\left(n+1\right)\left(n+2\right){\left(\frac{R}{r}\right)}^{n+1}{P}_n\left(\text{cos} \psi \right)}$$

(17)

A practical computation of the Moho density contrast is realized in two numerical steps. The unknown parameters W are first computed by solving the inverse to the system of observation equations. These observation equations are formed by applying a discretization to the integral equation in Equation (16). The estimated parameters W are then used to compute the Moho density contrast according to the definition given in Equation (12).

2.3. Combined Model

The functional model from Equation (16) is now used to formulate a method for finding the Moho density contrast by combining the gravity-gradient data and seismic crustal model. For this purpose we use the a priori information on the Moho depth and density contrast from an available seismic crustal model to define the condition equation for solving the VMM problem. A least-squares technique is then applied to estimate the Moho density contrast from the gravity-gradient data. Since the accuracy of seismic data is typically not provided, we do not apply a stochastic model.

The condition equation is defined in the following form

$$\Delta \rho T-\Delta {\rho }_0{T}_0+M=0$$

(18)

where the parameter M is defined by

$$M=\frac{1}{4\pi G}\left[{\displaystyle \sum _{n=0}^{\infty }\frac{2n+1}{n+1}}\left(-g_n^{\mathrm{T}}-g_n^{\mathrm{B}}-g_n^{\mathrm{S}}\right)+W\right]$$

(19)

The condition equation in Equation (18) assumes that the differences (i.e., residuals) between the gravimetric and seismic Moho parameters are minimized by means of applying a least-squares technique. The system of condition equations then becomes

$$\mathbf{B}=\left[\begin{array}{ccccc}T& \Delta \rho & -T_0& -\Delta {\rho }_0& 1\end{array}\right]$$

(20)

and

$$\mathbf{w}=0-\Delta \rho T+\Delta {\rho }_0{T}_0+M$$

(21)

The system of condition equations can uniquely be solved by applying the minimum-norm condition. Hence, we have

$$\widehat{v}=\left(\begin{array}{c}\delta \widehat{T}\\ \delta \Delta \widehat{\rho }\\ \delta {\widehat{T}}_0\\ \delta \Delta {\widehat{\rho }}_0\\ \delta \widehat{M}\end{array}\right)={\mathbf{B}}^{\mathrm{T}}{\left(\mathbf{B} {\mathbf{B}}^{\mathrm{T}}\right)}^{-1}\mathbf{w}$$

(22)

The solution of Equation (22) yields the correction terms to the Moho parameters.

2.4. Kernel Behavior

To compute the parameter W, the integral equation in Equation (16) is discretized and then solved numerically. The dependence of result on the integration domain can be investigated from a spatial behavior of the integral kernel L. The gravity inversion could be localized if the integral kernel attenuates quickly with an increasing spherical distance [34,35]. The inversion is then carried out using only data within the near zone, while the far-zone contribution is disregarded. The near zone is chosen so that truncation errors due to disregarding the far-zone contribution are negligible.

A spatial behavior of the integral kernel L is illustrated in Figure 1. The kernel attenuates very quickly already at 2 arc-deg of the spherical distance with nearly zero values at 4 arc-deg. The near-zone limit ψ ≤ 4 arc-deg thus sufficiently reduces truncation errors.

3. Numerical Realization

In this section we begin with a brief tectonic classification of the study area, then specify input data and applied methods and finally present and validate results.

3.1. Study Area and Tectonic Setting

The regional study area is bounded by the parallels 0 and 70 arc-deg of the northern latitude and the meridians 40 and 120 arc-deg of the eastern longitude. A tectonic configuration of this study area comprises parts of the Eurasian, Indian, Arabian, African, Amur and Sunda plates including the Tibetan, Iranian and Yangtze tectonic blocks (Figure 2a). It comprises zones of the compressional, extensional and strike-slip tectonism. The most prominent tectonic feature is the active continent-to-continent collision of the Indian plate with the Tibetan block, which has been responsible for the uplift of Himalaya. Tomographic evidence also suggests the southward subduction of the Eurasian lithosphere beneath the Tibetan block. A significant compressional tectonism due to the subduction of the Indian and Eurasian lithosphere beneath the Tibetan block resulted in the uplift of the whole Tibetan Plateau. The continent-to-continent collision between the Arabian and Eurasian plates uplifted the Iranian block, which is separated from the Eurasian plate by the strike-slip fault systems with major sections of the Herat, Chakaneh-Neyshabur, Bakharden-Quchan, Kopet Dagh and Mosha faults. The extensional tectonism occurred along the Mid-Indian oceanic rift zone and formed also the Baikal Rift Zone. Another significant feature within the study area is the Ural Orogen, which represents the oldest known orogenic structure mostly deeply buried into younger sediments. The major geological provinces include also continental basins (Indo-Ganges Basin, Tarim Basin and Sichuan Basin) and cratons (North China Carton, Yangtze Craton, Siberian Caron and Yakutai Craton).

3.2. Data Acquisition

The topographic heights on land and the bathymetric depths offshore from the SRTM30 database were used to compute the topographic and bathymetric-stripping gravity corrections with a spectral resolution complete to the spherical harmonic degree 180. These gravity corrections were computed on a 1 × 1 arc-deg surface grid. The average density of the upper continental crust 2670 kg·m⁻³ [36] was used to compute the topographic gravity correction. The bathymetric-stripping gravity correction was evaluated by applying the seawater density-depth equation [37]. The sediment-stripping gravity correction was calculated using the CRUST1.0 sediment data. The solid topography of the study area is shown in Figure 2a. The combined topographic-bathymetric-sediment gravitational contribution within the study area varies between −354 and 753 mGal (see Figure 2b). This contribution is mostly negative offshore, while positive on land with maxima over Himalaya and Tibet.

We applied the combined method (Section 2.3) to estimate the Moho density contrast from the GOCE Level 2 EGG_TRF_2 data product [38] over a period from 1st of January until 31st of May 2012 and the CRUST1.0 seismic crustal model. The GOCE data defined in local north oriented frame (LNOF) were transformed into the vertical gravity-gradient gradients (in the geocentric reference frame) and reduced for the normal vertical gravity-gradient computed according to the GRS80 parameters [39]. The GOCE vertical gravity gradients (at respective satellite altitudes) reach positive values up to 1.3 Eötvös (E) as well as negative values (−1.4 E), with the largest horizontal spatial variations across profiles intersecting Tarim and Ganges Basins and Tibetan Orogen (see Figure 3).

3.3. Estimation of the Parameter W

The GOCE vertical gravity gradients V_rr (Figure 3) were first inverted into the parameters W according to Equation (16). For this purpose, the integral equation was discretized and the (unknown) parameters W were then computed by applying the (least-squares) spherical harmonic analysis. The system of observation equations was solved directly by a matrix inversion and the (zero-order) Tikhonov regularization was applied to stabilize the ill-posed problem. The regularization parameter was estimated by applying the quasi-optimal method [40]. The computed parameters W vary between −21 and 24 mGal (see Figure 4). Since the values V_rr and W are highly correlated, their spatial patterns are very similar (see Figure 3 and Figure 4).

3.4. Combined Solution for the Moho Parameters

We further computed the Moho depth and density contrast within the study area on a 1 × 1 arc-deg spherical grid by solving the system of condition equations in Equation (22). This solution combines the parameters W estimated from the vertical gravity gradients, the combined topographic-bathymetric gravitational contribution computed from the SRTM30 data, the gravitational contribution of sediments estimated using the CRUST1.0 sediment data and the CRUST1.0 Moho parameters T₀ and ∆ρ₀. For the least-squares estimation, the weight matrix was chosen to be the identity matrix. We then repeated the computation for the Moho parameter T₀ taken from the seismic models GSMM [41] and M13 [42]. Since these two models provide information only on the Moho depth, we adopted the value ∆ρ₀ from CRUST1.0. The combined solutions (denoted as: R1_CRUST1.0, R2_GSMM and R3_M13) are shown in Figure 5 and their statistical summaries are given in

Table 1. Statistics of the Moho depth solutions.

Combined Solution	Min (km)	Max (km)	Mean (km)	STD (km)
R1_CRUST1.0	7.4	62.4	33.2	10.3
R2_GSMM	7.7	64.2	32.2	9.9
R3_M13	10.2	61.6	34.5	10.3

(for the Moho depth) and

Table 2. Statistics of the Moho density contrast solutions.

Combined Solution	Min (kg·m−3)	Max (kg·m−3)	Mean (kg·m−3)	STD (kg·m−3)
R1_CRUST1.0	166	624	166	87
R2_GSMM	211	602	211	75
R3_M13	149	640	149	87

(for the Moho density contrast).

As seen in Figure 5, all three results exhibit a similar spatial pattern in estimated values of the Moho depth and density contrast. The most pronounced feature in the Moho geometry is the contrast between a thin oceanic crust compared to a thicker continental crustal structure with the maximum Moho deepening under orogens of Himalaya and Tibet. The most pronounced pattern in the Moho density contrast is related with significant differences between the density contrast under continental and oceanic crustal structures. In this case, however, the density contrast under orogens is not significantly enhanced. Instead, the Moho density contrast along mid-oceanic rift zones is more pronounced with respect to the density contrast under the remaining oceanic crust. As seen, the Moho density contrast under the oceanic crust varies typically from 200 to 400 kg·m⁻³ with minima along mid-oceanic rift zones. The Moho density contrast under the continental crust mostly exceeds 450 kg·m⁻³. Despite these similarities between all three combined results, some more localized differences between them could be recognized. However, these differences are within the level of expected uncertainties of the estimated Moho parameters. It was demonstrated, for instance, by [43] that the Moho depth uncertainties (estimated using seismic data) beneath Europe could regionally exceed 10 km, with an average error about 4 km. Even larger errors are to be considered at regions where seismic data are sparse or absent. Similarly, relative errors 10% or larger could be expected in estimated values of the Moho density contrast.

In our numerical realization we disregarded the density heterogeneities within the remaining crust down to the Moho interface, because the CRUST1.0 data of consolidated crustal density layers are not sufficiently accurate. To demonstrate this, we computed the combined gravitational contribution of topography, bathymetry, sediments and consolidated crust. This contribution varies between −514 to 943 mGal, with a mean 198 mGal and a standard deviation 265 mGal (see Figure 6a). We then used these values to estimate the Moho depth and density contrast. The application of the consolidated crust-stripping gravity correction introduced a systematic bias in estimated values of the Moho depth. As seen in Figure 6b, the Moho depth in this case varies between 1.6 and 71.2 km, with a mean 41.3 km and a standard deviation 11.1 km. This mean Moho depth differs significantly from the Moho solutions obtained without applying the consolidated crust-stripping gravity correction (cf. Table 1). The Moho density contrast is shown in Figure 6c. It varies between 159 and 712 kg·m⁻³, with a mean 537 kg·m⁻³ and a standard deviation 65 kg·m⁻³. The comparison of these values with the Moho density contrast estimated without applying the consolidated crust-stripping gravity correction (cf. Table 2) indicates again the presence of a large systematic bias.

3.5. Validation of Results

To validate the combined results presented in Figure 5, we compared the Moho density contrast with the CRUST1.0 and KTH1.0 [44] models (see Figure 7 and statistics in

Table 3. Statistics of the Moho density contrast from the CRUST1.0 and KTH1.0 models.

Model	Min (kg·m−3)	Max (kg·m−3)	Mean (kg·m−3)	STD (kg·m−3)
CRUST1.0	88	694	495	95
KTH1.0	−40	610	374	94

). In addition, we compared the Moho depth with the CRUST1.0, KTH1.0, GSMM, M13 and GEMMA [45] models (see Figure 8 and statistics in

Table 4. Statistics of the Moho depth from the CRUST1.0, KTH1.0, GSMM, M13 and GEMMA models.

Model	Min (km)	Max (km)	Mean (km)	STD (km)
CRUST1.0	7.7	72.6	34.1	12.8
KTH1.0	6.2	85.8	35.1	13.8
GSMM	7.0	78.9	32.6	11.4
M13	8.7	74.4	35.4	13.2
GEMMA	5.5	69.7	33.5	10.1

). The Moho density contrast differences are plotted in Figure 9 and their statistical summary is given in

Table 5. Statistics of the Moho density contrast differences.

Differences	Min (kg·m−3)	Max (kg·m−3)	Mean (kg·m−3)	RMS (kg·m−3)
R1_CRUST1.0 − CRUST1.0	−253	200	−13	43
R1_CRUST1.0 − KTH1.0	−112	292	106	119
R2_GSMM − CRUST1.0	−253	200	−13	42
R2_GSMM − KTH1.0	−98	301	106	118
R3_M13 − CRUST1.0	−227	245	−11	56
R3_M13 − KTH1.0	−157	353	109	125

. For completeness, we also provided a statistical summary for the Moho depth differences in

Table 6. Statistics of the Moho depth differences.

Differences	Min (km)	Max (km)	Mean (km)	RMS (km)
R1_CRUST1.0 − CRUST1.0	−20.1	23.1	−0.1	5.5
R1_CRUST1.0 − KTH1.0	−25.2	19.8	−1.2	6.4
R1_CRUST1.0 − GSMM	−21.8	25.8	1.2	5.7
R1_CRUST1.0 − M13	−23.9	22.5	−1.5	6.0
R1_CRUST1.0 − GEMMA	−23.3	13.0	0.4	3.9
R2_GSMM − CRUST1.0	−23.3	21.5	−1.9	5.9
R2_GSMM − KTH1.0	−35.5	18.3	−2.8	6.9
R2_GSMM − GSMM	−29.7	18.3	−0.3	5.5
R2_GSMM − M13	−29.1	21.2	−3.1	6.7
R2_GSMM − GEMMA	−26.3	11.4	−1.3	4.3
R3_M13 − CRUST1.0	−18.8	24.3	1.1	5.7
R3_M13 − KTH1.0	−27.9	21.1	0.1	6.3
R3_M13 − GSMM	−21.5	25.9	2.5	6.2
R3_M13 − M13	−21.0	23.9	−0.1	5.7
R3_M13 − GEMMA	−22.0	16.2	1.7	4.3

.

As seen in Table 5, all three combined solutions for the Moho density contrast better agree (by means of the RMS of differences) with the CRUST1.0 seismic model than with the KTH1.0 combined gravimetric-seismic model. Moreover, the R2_GSMM combined solution relatively poorly reproduces the Moho density contrast under Himalaya and most of Tibet (cf. Figure 9b). These large density contrast differences are explained by an unrealistically shallow Moho of the GSMM seismic model (see Figure 8c) under these orogens. Despite the R2_GSMM combined solution has the best RMS fit with the CRUST1.0 and KTH1.0 models (cf. Table 5), this combined solution is likely very inaccurate in these areas. The R1_CRUST1.0 and R3_M13 combined models are, on the other hand, very similar by means of their spatial pattern as well as their RMS fit with testing models. This finding is also confirmed by the results of the Moho depth modeling. As seen in Table 6, the RMS fit of these two combined solutions with the Moho depth from the CRUST1.0, KTH1.0, M13 and GEMMA models are similar. Interestingly, however, the GSMM Moho depth errors in Himalaya and Tibet are to a large extent reduced in the R2_GSMM combined solution. This is due to the fact that the condition equations in Equation (18) were specified for the mean values of the Moho parameters T₀ and ∆ρ₀. The formulation of condition equations for mean Moho parameters has also practical implications, because in this way we could to some extent reduce errors in estimated Moho parameters over areas with insufficient seismic-data coverage.

4. Summary and Concluding Remarks

We have developed the combined method for a determination of the Moho parameters from the GOCE in-orbit gravity gradients and constrained using seismic crustal model. We applied this method to investigate the Moho density contrast under most of Eurasia including surrounding continental and oceanic areas. The input data used for computation involved the topographic information from the SRTM on land, the bathymetric information derived from processing the satellite-altimetry data offshore, the sediment density and thickness information from the CRUST1.0 data and the vertical gravity gradients obtained from processing the GOCE satellite observables. The least-squares technique, applied for a regional Moho recovery, was formulated for the system of condition equations and solved based on finding the minimum-norm solution in order to obtain an optimal fit between gravimetric and seismic models. The condition equations were defined for mean values of seismic Moho parameters in order to suppress errors from areas with insufficient seismic-data coverage.

Our results confirmed large variations of the Moho density contrast with minima along the mid-oceanic rift zones and maxima under most of continental crustal structures. These results have a relatively good agreement with the CRUST1.0 and KTH1.0 models, but differ substantially from some previously published results of gravimetric studies. The most significant inconsistency was found in overestimated values of the Moho density contrast under Tibet and Himalaya, which according to [23,30] could reach 800 kg·m⁻³ or even more. In more recent study [43], however, much smaller values of the Moho density contrast were estimated (i.e., the KTH1.0 model). Moreover, the Moho density contrast according to [46] could hardly exceed 600 kg·m⁻³ and much smaller values of the Moho density contrast were also reported in various seismic studies. In this study, we have confirmed these findings by showing that the Moho density contrast under most of the Eurasian continental crust is typically very uniform with most of the values being within a range from 400 to 600 kg·m⁻³.