Moho Density Contrast in Central Eurasia from GOCE Gravity Gradients

Mehdi Eshagh 1, Matloob Hussain 1,2, Robert Tenzer 3,4,* and Mohsen Romeshkani 5 1 Department of Engineering Science, University West, Trollhättan 46186, Sweden; mehdi.eshagh@hv.se (M.E.); mutloob.hussain@hv.se (M.H.) 2 Department of Earth Sciences, Quaid-i-Azam University, Islamabad 45320, Pakistan 3 The Key Laboratory of Geospace Environment and Geodesy, Wuhan University, Wuhan 430079, China 4 New Technologies for the Information Society (NTIS), University of West Bohemia, Plzen 30614, Czech Republic 5 School of Surveying and Geospatial Engineering, College of Engineering, University of Tehran, Tehran 14395-515, Iran; romeshkani@ut.ac.ir * Correspondence: rtenzer@sgg.whu.edu.cn; Tel.: +86-27-6877-8649


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 ´3 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 ´3 (for the mafic lower crust) to 440 kg¨m ´3 (for the felsic lower crust), with an apparently typical value 440 kg¨m ´3 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 ´3, while this range is between 340 and 790 kg¨m ´3 when computed with respect to the reference crustal density 2670 kg¨m ´3 (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 ´3 (in the Pacific region) to 988 kg¨m ´3 (beneath the Tibetan Plateau).A similar range of values between 82 and 965 kg¨m ´3 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.0values 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.

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.

Moho Density Contrast from Gravity Disturbance
The VMM isostatic problem was formulated in the following generic form [26,27,32] where G is the Newton's gravitational constant, R is the Earth's mean radius, ∆ρ is the (variable) Moho density contrast, T is the Moho depth, δg i is the isostatic gravity disturbance, σ is the unit sphere, and dσ 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 pψ, Tq " where P n is the Legendre polynomial of degree n.The isostatic gravity disturbances δg i on the right-hand side of Equation ( 1) are computed from the gravity disturbances δg by subtracting the gravitational contributions of topography g T , bathymetry g B and sediments g S and adding the compensation attraction A C0 .Hence The compensation attraction A C0 in Equation ( 3) is computed from [26] A C0 " where T 0 and ∆ρ 0 are mean values of the Moho depth and density contrast respectively.To solve the VMM problem for finding the Moho density contrast ∆ρ, we first simplified the integral kernel K in Equation (2).By applying a Taylor series (up to a first-order term) to p1 ´T{Rq n`3 , we get K pψ, Tq " T R The surface integral on the left-hand side of Equation ( 6) is defined in terms of the Laplace harmonic p∆ρTq n as follows [33] σ p∆ρTq P n pcos ψq dσ " 4π 2n `1 p∆ρTq n Inserting from Equation (7) to Equation (6), we obtained the following relation The coefficients p∆ρTq n in Equation ( 8) are generated from the Laplace harmonics δg i n of the isostatic gravity field as follows where g T n , g B n and g S n are, respectively, the Laplace harmonics of the topographic, bathymetric and sediment gravitational contributions, and δ n0 is Kronecker's delta.
The expression in Equation ( 9) is finally rearranged into the following form 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.

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 δg n reads 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 where the parameter W is defined by 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 Alternatively, the Laplace harmonics W n can be defined in the following form [33] W Substituting from Equation ( 15) back to Equation ( 14) and applying further simplifications, we arrived at R 4π σ W L pr, ψq dσ " r 2 V rr (16) where the integral kernel L reads L pr, ψq " 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).

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 where the parameter M is defined by 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 and The system of condition equations can uniquely be solved by applying the minimum-norm condition.Hence, we have The solution of Equation ( 22) yields the correction terms to the Moho parameters.

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.

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.

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-tocontinent 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).

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 −3 [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.0sediment data.The solid topography of the study area is shown in Figure 2a.The combined topographic-bathymetricsediment gravitational contribution within the study area varies between −354 and 753 mGal (see

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.

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).

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 ´3 [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.0sediment 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.0seismic 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).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.0seismic 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).

Estimation of the Parameter W
The GOCE vertical gravity gradients Vrr (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 Vrr and W are highly correlated, their spatial patterns are very similar (see Figures 3 and 4).

Combined Solution for the Moho Parameters
We further computed the Moho depth and density contrast within the study area on a 1 × 1 arcdeg 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-

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 Figures 3 and 4).

Estimation of the Parameter W
The GOCE vertical gravity gradients Vrr (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 Vrr and W are highly correlated, their spatial patterns are very similar (see Figures 3 and 4).

Combined Solution for the Moho Parameters
We further computed the Moho depth and density contrast within the study area on a 1 × 1 arcdeg 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-

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.0sediment data and the CRUST1.0Moho parameters T 0 and ∆ρ 0 .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 0 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 ∆ρ 0 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 (for the Moho depth) and bathymetric gravitational contribution computed from the SRTM30 data, the gravitational contribution of sediments estimated using the CRUST1.0sediment data and the CRUST1.0Moho parameters T0 and ∆ρ0.For the least-squares estimation, the weight matrix was chosen to be the identity matrix.We then repeated the computation for the Moho parameter T0 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 ∆ρ0 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 (for the Moho depth) and Table 2 (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 ´3 with minima along mid-oceanic rift zones.The Moho density contrast under the continental crust mostly exceeds 450 kg¨m ´3.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.0data 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 ´3, with a mean 537 kg¨m ´3 and a standard deviation 65 kg¨m ´3.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.

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).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).The Moho density contrast differences are plotted in Figure 9 and their statistical summary is given in Table 5.For completeness, we also provided a statistical summary for the Moho depth differences in Table 6.

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).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).The Moho density contrast differences are plotted in Figure 9 and their statistical summary is given in Table 5.For completeness, we also provided a statistical summary for the Moho depth differences in Table 6.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.0seismic 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 T0 and ∆ρ0.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.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.0seismic 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 0 and ∆ρ 0 .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.

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.0data 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 ´3 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 ´3 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 ´3.

A 19 Figure 1 .
Figure 1.Integral kernel L (for the mean satellite attitude 250 km).

Figure 1 .
Figure 1.Integral kernel L (for the mean satellite attitude 250 km).

FigureFigure 2 .
Figure 2b).This contribution is mostly negative offshore, while positive on land with maxima over Himalaya and Tibet.

Figure 2 .
Figure 2. Regional maps of: (a) the solid topography and (b) the combined topographic-bathymetric gravitational contribution.Tectonic plate boundaries are indicated by red line.

Figure 4 .
Figure 4.Estimated values of the parameter W.

Figure 4 .
Figure 4.Estimated values of the parameter W.

Figure 4 .
Figure 4.Estimated values of the parameter W.

Figure 6 .
Figure 6.Effect of the consolidated crust-stripping gravity correction to Moho parameters: (a) the combined gravitational contribution of topography, bathymetry, sediments and consolidated crust.The estimated values of: (b) Moho depth and (c) Moho density contrast.

Figure 6 .
Figure 6.Effect of the consolidated crust-stripping gravity correction to Moho parameters: (a) the combined gravitational contribution of topography, bathymetry, sediments and consolidated crust.The estimated values of: (b) Moho depth and (c) Moho density contrast.

Table 2 (
for the Moho density contrast).

Table 1 .
Statistics of the Moho depth solutions.

Table 2 .
Statistics of the Moho density contrast solutions.

Table 5 .
Statistics of the Moho density contrast differences.

Table 6 .
Statistics of the Moho depth differences.