Reformulation of the Vening-Meinesz Moritz Inverse Problem of Isostasy for Isostatic Gravity Disturbances

The isostatic gravity anomalies have been traditionally used to solve the inverse problems of isostasy. Since gravity measurements are nowadays carried out together with GPS positioning, the utilization of gravity disturbances in various regional gravimetric applications becomes possible. In global studies, the gravity disturbances can be computed using global geopotential models which are currently available to a relatively high accuracy and resolution. In this study we facilitate the definition of the isostatic gravity disturbances in the Vening-Meinesz Moritz inverse problem of isostasy for finding the Moho depths. We further utilize uniform mathematical formalism in the gravimetric forward modelling based on methods for a spherical harmonic analysis and synthesis of gravity field. We then apply both mathematical procedures to determine globally the Moho depths using the isostatic gravity disturbances. The results of gravimetric inversion are finally compared with the global crustal seismic model CRUST2.0; the RMS fit of the gravimetric Moho model with CRUST2.0 is 5.3 km. This is considerably better than the RMS fit of 7.0 km obtained after using the isostatic gravity anomalies.


Introduction
The functional models of solving the inverse problem of isostasy have been traditionally formulated by means of the isostatic gravity anomalies (cf.Heiskanen and Moritz [1], p. 133).According to these definitions the topographic mass surplus and the ocean mass deficiency are compensated either by a variable crustal thickness or density.The isostatic equilibrium is described in terms of the isostatic gravity anomalies which should theoretically be equal zero, provided that the refined Bouguer gravity anomalies are isostatically fully rebalanced by the corresponding gravitational attraction of compensating masses.The Pratt-Hayford isostatic model is based on adopting a constant depth of compensation while considering a variable density contrast (Pratt [2], Hayford [3], Hayford and Bowie [4]).In the Airy-Heiskanen isostatic model a constant density contrast is assumed while a depth of compensation is variable (Airy [5], Heiskanen and Vening Meinesz [6]).Both these classical isostatic models are based on a local compensation scheme.Vening-Meinesz [7] modified the Airy-Heiskanen theory of isostasy with applying the regional instead of local compensation.Moritz [8] utilized the Vening-Meinesz inverse problem of isostasy for finding the Moho depths.Sjöberg [9] fur-ther generalized this concept for finding the Moho depths and density contrast.Sjöberg and Bagherbandi [10] computed the Moho depths based on solving Moritz's generalization of the Vening-Meinesz inverse problem of isostasy (VMM isostatic model).Bagherbandi and Sjöberg [11] demonstrated that the VMM Moho depths better agree with the Moho data taken from the global crustal seismic model CRUST2.0(Bassin et al. [12]) than those obtained based on solving the Airy-Heiskanen isostatic model.
The results of regional and global studies have shown often existing significant disagreement between the isostatic and seismic Moho depths.Several different reasons explaining this misfit were proposed and also confirmed by numerical experiments.Kaban et al. [13], for instance, demonstrated that the isostatic compensation does not take place only within the Earth's crust but essentially also within the lithospheric mantle.This finding was later confirmed by Kaban et al. [14] and Tenzer et al. [15,16].Several authors argued that the isostatic balance is also partially affected by the changing rigidity, glacial isostatic adjustment, plate motion, ocean lithosphere cooling, and other geophysical processes.Moreover, large portion of the isostatic mass balance is attributed to variable li-thospheric density structures which are usually not taken into consideration in computing the isostatic gravity anomalies.Therefore, the models for gravimetric recovery of the Moho parameters should incorporate all known information on subsurface density structures.One example can be given in Greenland and Antarctica where the application of the ice density contrast stripping correction to gravity data is essential for a realistic interpretation of gravimetric results.Another significant gravitational contribution to be modelled and subsequently subtracted from gravity data is due to large sedimentary basins.Braitenberg et al. [17] and Wienecke et al. [18] demonstrated that the misfit of the isostatic assumption of the Moho interface to the long-wavelength part of gravity field is explained by large sedimentary basins and rigidity variations of the crustal plate.
The gravity anomalies have been primarily used in regional and global studies investigating the Earth's inner structures.Vajda et al. [19], however, argued that the definition of gravity disturbances in the context of these studies is theoretically more appropriate.Moreover, modern gravity observation techniques (such as air-born gravimetry) nowadays incorporate the GPS positioning systems.Therefore, it is expected that the gravity disturbances will become the most often used gravity data type in all gravimetric applications.This is due to the fact that GPS observations provide the geodetic heights above the reference ellipsoid surface, while the definition of gravity anomalies requires topographic heights with respect to sea level.Tenzer et al. [20][21][22] utilized the definition of gravity disturbances in the forward modeling of gravitational field generated by major known crustal density structures.Following this concept, here we define the VMM inverse problem of isostasy by means of the isostatic gravity disturbances.

Refined Gravity Disturbances
To solve the gravimetric problem of isostasy for finding the Moho parameters, the gravitational contributions of all known density contrasts within the Earth's crust should be modeled and subsequently removed from gravity data.Moreover, the inhomogeneous density structures within the mantle lithosphere and sub-lithosphere mantle should be taken into consideration provided that reliable data of global mantle density structures are available.The resulting gravity data which have a maximum correlation with the Moho geometry are theoretically the most appropriate for the gravimetric recovery of the Moho parameters.Tenzer et al. [23] used the global gravity and crustal density structure models to investigate the global correlation of various gravity field quantities with the Moho geometry.They demonstrated that a maximum correlation is attained when using the gravity distur-bances which are corrected for the gravitational contributions of topography and density contrasts of ocean, ice and sediments.They also showed that the application of the stripping gravity correction due to the anomalous density structures within the consolidated (crystalline) crust slightly decreased the correlation with the Moho geometry.The possible reason is more likely due to large inaccuracies within the CRUST2.0consolidated crustal data.cs g The gravity disturbances  which have theoretically a maximum correlation with the Moho geometry (for a chosen lithosphere density model) are obtained from the gravity disturbances g  after subtracting the gravitational effect of topographic masses.Furthermore, the gravitational contributions of density differences between the known and synthetic model of lithosphere are modeled and subsequently removed from the topography-corrected gravity disturbances.This is done by means of applying the stripping gravity corrections.The computation of the gravity disturbances and gravity corrections is done here in spectral domain using methods for a spherical harmonic analysis and synthesis of gravity field.

 
, r g  at a point The gravity disturbance  is computed using the following expression [1]  is the maximum degree of spherical harmonics.The coefficients , n m are obtained from the global geopotential model (GGM) coefficients after subtracting the spherical harmonic coefficients describing the normal gravity field.The 3-D position is defined in the system of spherical coordinates ; where is the spherical radius, and denotes the spherical direction with the spherical latitude  and longitude  .
Tenzer et al. [22] developed and applied the uniform mathematical formalism of computing the topographic and density contrasts stripping gravity corrections.It utilizes the expression for computing the gravitational attraction g (defined approximately as a negative radial derivative of the respective potential ; i.e., V g V r     ) generated by an arbitrary volumetric mass layer with a variable depth and thickness while having a laterally distributed vertical density variation.According to their approach, the gravity corrections at a point The potential coefficients in Equation (2) read where Earth 5500 kg/m 3 is the adopted value of the Earth's mean mass density.The numerical coefficients are defined as follows Fl , Fu : Fl The terms and define the spherical lower-and upper-bound laterally distributed radial density variation functions L and of degree .These spherical functions and their higher-order terms   and The infinitesimal surface element on the unit sphere is denoted as For a specific volumetric layer, the mass density  is either constant  , laterally-varying or-in the most general case-approximated by the laterally distributed radial density variation model using the following polynomial function (for each lateral column) where U is the nominal value of the lateral density stipulated at the depth U of the upper bound of the volumetric mass layer.This density distribution model describes the radial density variation within the volumetric mass layer at a location .Alternatively, when modeling the gravitational field of the anomalous mass density structures, the density contrast     of the volumetric mass layer relative to the reference crustal density c  is defined as where U is the nominal value of the lateral density contrast stipulated at the depth of the upper bound of the volumetric mass layer.
The coefficients , n m and , U n m combine the information on the geometry and density (or density contrast) distribution of volumetric layer.These coefficients are generated to a certain degree of spherical harmonics using the discrete data of the spatial density distribution (i.e., typically provided by means of density, depth and thickness data) of a particular structural component of the Earth's interior.the constant density values of the Earth's crust and the encompassing upper(most) mantle respectively (cf.Vening Meinesz [7]).The isostatic gravity disturbance

Vening Meinesz-Moritz Isostatic Model
where g  is the refined gravity disturbance (defined in Section 2), and c g is the gravitational attraction of isostatic compensation masses (e.g., Moritz [8]).
Sjöberg [9] derived the VMM inverse problem of isostasy in the following generic form where m 3 •kg -1 •s -2 is Newton's gravitational constant.The integral kernel K in Equation ( 11) is a function of parameters G  and s ; where  is the spherical distance between two points   , r  and , and where n is the Legendre polynomial of degree .The isostatic gravity disturbance functional f on the righthand side of Equation ( 11) is introduced as follows The expression given in Equation ( 11) is a non-linear Fredholm integral equation of the first kind.A direct solution for finding the Moho depths was given by Sjöberg [9] which is a second-order approximation.It reads The term 1 is computed in spectral domain using the following expression The numerical coefficients , n m f of the isostatic gravity disturbance functional f are defined as stipulated at the sphere of radius is given by (cf.Sjöberg [9]) R where 0 is the adopted nominal mean value of the Moho depths.

Input Data Acquisition
The expressions for the gravimetric forward modeling were utilized to compute the topographic and stripping gravity corrections due to ocean, ice and sediment density contrasts.These corrections were then applied to the gravity disturbances in order to obtain the refined gravity disturbances which were used for solving the VMM inverse problem of isostasy.All computations were realized globally on a 1 × 1 arc-deg geographical grid of the surface points.The gravity disturbances were generated using the EGM08 coefficients (Pavlis et al. [24]) with a spectral resolution complete to degree 180 of spherical harmonics (which corresponds to a half-wavelength of 1 arc-deg or about 100 km).The spherical harmonic terms of the normal gravitational field were computed according to the parameters of GRS-80 (Moritz [25]).The same spectral resolution was used to compute the topographic and bathymetric (ocean density contrast) stripping gravity corrections.These two gravity corrections were computed from the DTM2006.0 coefficients (Pavlis et al. [26]).The global topographic/bathymetric model DTM-2006.0was released together with EGM2008 by the US National Geospatial-Intelligence Agency EGM development team.The average density of the upper continental crust 2670 kg/m 3 (cf.Hinze [27]) was adopted for defining the topographic and reference crustal densities.The bathymetric stripping gravity correction was computed using the depth-dependent seawater density model (see Tenzer et al. [28]).This empirical ocean density model was developed by Gladkikh and Tenzer [29] based on the analysis of oceanographic data of the World Ocean Atlas 2009 (WOA09) and the World Ocean Circulation Experiment 2004 (WOCE04).WOA09 oceanographic products are made available by NOAA's National Oceanographic Data Center (Johnson et al. [30]).The WOCE04 datasets are provided by the German Federal Maritime and Hydrographic Agency (Gouretski and Koltermann [31]).Tenzer et al. [32] acquired, based on the comparison of the experimental and theoretical seawater density values, that this empirical model approximates the actual seawater density distribution with the maximum relative error better than 0.6%, while the corresponding average error is about 0.1%.For the adopted value of the reference crustal density 2670 kg/m 3 and surface seawater density 1027.91 kg/m 3 (cf.Gladkikh and Tenzer [29]) the nominal ocean density contrast (at zero depth) equals 1642.09kg/m 3 .The parameters of the depth-dependent ocean density term in Equation ( 9    contrast stripping gravity correction was computed individually for the CRUST2.0soft and hard sediments.The results are shown in Figures 5(a) and (b).The soft sediment stripping gravity correction globally varies within 7 and 144 mGal, with a mean of 33 mGal and a standard deviation is 21 mGal.The corresponding correction due to the hard sediment density contrast is within -7 and 89 mGal, with a mean of 11 mGal and a standard deviation is 13 mGal.The complete sediment density contrast stripping gravity correction is everywhere positive and varies within 14 and 125 mGal, with a mean of 35 mGal and a standard deviation is 20 mGal.
The refined gravity disturbances obtained after applying the topographic and stripping gravity corrections due to the ocean, ice and sediment density contrasts are shown in Figure 6.These refined gravity data globally vary between −498 and 760 mGal, with a mean of 320 mGal and a standard deviation is 196 mGal.Whereas the gravity disturbances globally vary mostly within ±300 mGal, the application of topographic and stripping gravity corrections changed the gravity field significantly.The range of refined gravity disturbances (of 1258 mGal) is more than twice as large as the range of uncorrected gravity disturbances (of 596 mGal).The largest gravity changes over continents are due to applying the topographic gravity correction especially in mountainous regions.The application of the bathymetric striping gravity correction changed substantially the gravity field over oceans.The ice density contrast stripping gravity correction changed the gravity field in central Greenland and Antarctica.Less substantial changes in gravity field were found after applying the sediment density contrast stripping gravity correction.The maxima of these changes are along the continental margins and in polar areas with the largest sediment deposits.The global gravity field obtained after step-wise application of these corrections  were shown in Tenzer et al. [35,36].

Isostatic Moho Recovery
The refined gravity data (which have a maximum correlation with the Moho geometry) shown in Figure 6 were further utilized in definition of the isostatic gravity disturbances (in Equation ( 10)).The gravitational contribution of crustal compensation masses was computed approximately using the expression in Equation ( 17) for the nominal compensation attraction.The isostatic gravity disturbances were then used to determinate the Moho depths on a 1 × 1 arc-deg global grid.The computation was carried out based on solving the VMM inverse problem of isostasy (see Section 3).The result is shown in Figure 7.The Moho depths globally vary between 2.1 and 61.5 km, with a mean of 22.9 km and a standard deviation is 12.1 km.The Moho depths computed and presented in Section 5 were compared with CRUST2.0 seismic data.The differences between the VMM and CRUST2.0Moho depths are within -25.8 and 26.4 km, the mean and RMS of these differences are 0.1 and 5.3 km respectively.We further repeated the whole computation using the isostatic gravity anomalies.The comparison of this result, not shown herein, with CRUST2.0Moho depths showed that the mean and RMS of differences between them are -3.4 and 7.0 km respectively.The range of differences is within -26.6 and 28.0 km.The result based on using the isostatic gravity disturbances thus better agrees with the CRUST2.0Moho depths by means of the RMS fit.A significant improvement was also achieved by minimizing the systematic bias between these solutions.
As seen from these results, relatively small differences between the gravity disturbances and gravity anomalies correspond to relatively large differences between the Moho results obtained based on using these two gravity data types.The differences between the EGM08 gravity disturbances and gravity anomalies computed with a spectral resolution complete to degree 180 of spherical harmonics are shown in Figure 8.These differences are within -33 and 26 mGal, the mean and RMS of these differences are -0.2 mGal and 9 mGal respectively.The corresponding differences between the Moho depths computed using these two gravity data types are shown in Figure 9.These differences are within -1.5 and 9.3 km, the mean and RMS of these differences are 2.3 km and 3.3 km respectively.

Discussion
The refined gravity disturbances used as the input gravity data for finding the Moho depths should (optimally) comprise only the gravitational signal of the Moho geometry.The global (absolute) correlation between these refined gravity data and the Moho geometry was confirmed to be 0.98.The currently available global models of the Earth's gravity field, topography, ice and bathymetry have a relatively high resolution and accuracy.The computation of the gravity data which are corrected for the gravitational contributions of these density components can thus be done with a sufficient accuracy.Large inaccuracies are, however, to be expected due to uncertainties of the currently available global crustal models.The current datasets of the spatial density distribution of sediment and consolidated (crystalline) crust have a low accuracy as well as resolution.As consequence, the computation of respective gravity corrections and corrected gravity data is still restricted.Tenzer et al. [16] estimated that the relative errors can reach as much as  10% in the refined gravity data.Moreover, the gravitational signals of the mantle lithosphere and sub-lithospheric structure (including the core-mantle geometry) are still presented in these refined gravity data.In the absence of reliable global models of mantle structures these gravitational contributions might be treated in the spectral domain by subtracting a long-wavelength part of gravity spectra (to a certain degree of spherical harmonics).This can be done while assuming that the subtracted long-wavelength gravity contribution is attributed mainly to the mantle structure and the core-mantle boundary.However, the current knowledge the spatial mantle density structure is restricted by the lack of reliable global data.A possible way how to estimate the maximum degree of long-wavelength spherical harmonics which should be removed from the refined gravity field was given by Eckhardt [37].The principle of this procedure is based on finding the representative depth of gravity signal attributed to each spherical harmonic degree term.
The spherical harmonics which have the depth below a certain limit (chosen, for instance, as the maximum Moho depth) are then removed from the gravity field.Nonetheless, the complete subtraction of the mantle gravity signal using this procedure is still questionable due to the fact that there is hardly any unique spectral distinction between the long-wavelength gravity signal from the mantle and the expected higher-frequency signal of the Moho geometry.Tenzer et al. [16] demonstrated the presence of significant correlation (>0.6) between the mantle gravity signal and the Moho geometry at the medium gravity spectrum (between 60 and 90 of spherical harmonics degree terms).On the other hand, the gravity signal of the core-mantle boundary and deep mantle structures could almost completely be subtracted from the refined gravity data as it is mainly attributed to the long-wavelength part of gravity spectra.
As seen from the functional model in Section 3, the errors in computed gravity data propagate proportionally to the Moho depth errors.The expected relative uncertainties in gravity data thus cause the errors of about 10% in the estimated Moho depths.Čadek and Martinec [38] estimated the uncertainties of Moho depths in their global crust thickness model to be about 20% (5 km) for the oceanic crust and of about 10% (3 km) for the continental crust.The results of more recent seismic and gravity studies, however, revealed that these error estimates are too optimistic.Grad et al. [39] demonstrated that the Moho uncertainties (estimated based on processing the seismic data) under the Europe regionally exceed 10 km with the average error of more than 4 km.Much large Moho uncertainties are obviously expected over large parts of the world where seismic data are absent or insufficient.A significant contribution to these Moho uncertainties is expected to be explained by several geophysical phenomena which are not accounted for in the isostatic functional model.Since the geophysical processes which contribute to isostatic imbalance have a different regional character, the global model which accounts for all these contributions (such as the glacial isostatic adjustment or oceanic lithosphere cooling) is difficult to establish and solve numerically.A possible method how to overcome to some extent these practical limitations was proposed and applied by Bagherbandi and Sjöberg [40].They combine the gravity and seismic data in order to model and subsequently account for the differences between the isostatic and seismic Moho models.The application of this method is out of the scope of this study.

Summary and Conclusions
We have redefined the Vening-Meinesz Moritz inverse problem of isostasy for the isostatic gravity disturbances while adopting the double layer model for defining the Moho density interface.The definition of the isostatic gravity disturbances was based on the refined gravity disturbances which have a maximum correlation with the Moho geometry.With reference to results of the correlation analysis by Tenzer et al. [23], the gravity disturbances corrected for the gravitational contributions of topography and density contrasts due to ocean, ice and sediments hold this condition.The application of the additional stripping gravity correction accounting for anomalous density structures within the crystalline crust was not taken into consideration as the currently available global data of this crustal component are still unreliable.Moreover, the mantle density structures were not modeled too due to the same reason.The methods for a spherical harmonic analysis and synthesis were used for computing the refined gravity disturbances.A spectral representation was also used in definition of the observation equation for solving the VMM model.The developed computational schemes were then applied to com-pute the isostatic gravity disturbances and to determine the Moho depths.The numerical experiment was carried out globally on a 1 × 1 arc-deg grid.The results were validated using the CRUST2.0Moho depths.
The global map in Figure 6 revealed the signature of gravity signal of which pattern corresponds with a spatial geometry of the Moho interface.The negative gravity values are prevailing over continents, while oceanic areas are dominated by the positive gravity values.The gravity minima agree with locations of the maximum continental crustal thickness under Himalayas, Tibet Plateau and Andes.The corresponding gravity maxima are seen especially in central Pacific Ocean.The contrast between the oceanic and continental lithospheric plate boundaries is marked along continental margins by the absolute gravity minima.The signature of the mantle lithosphere structures is also presented especially along the midoceanic ridges.
The Moho depths computed using the isostatic gravity disturbances based on solving the VMM model have a better agreement with the CRUST2.0seismic model than those computed using the isostatic gravity anomalies.The RMS fit of the VMM Moho depths with CRUST2.0 for the isostatic gravity disturbances was found to be 5.3 km.This RMS fit is about 24% better than the corresponding RMS fit of 7.0 km obtained when using the isostatic gravity anomalies.The facilitation of the isostatic gravity disturbances in the VMM model improved significantly the systematic bias otherwise found when using the isostatic gravity anomalies.The mean of differences between the VMM and CRUST2.0Moho depths obtained after using the isostatic gravity disturbances and gravity anomalies was found to be 0.1 and -3.4 km respectively.The systematic bias was thus almost completely eliminated.This numerical improvement is quite remarkable as the differences between the gravity anomalies and disturbances are globally mostly within 30 mGal.A possible explanation for such improvement might be due to the fact that the differences between these two gravity data types have a long-wavelength character.The respective changes in Moho results are thus more substantial than those attributed to the changes at a high-frequency part of gravity spectra.
Here we adopt the principle of solving Moritz's generalization of the Vening-Meinesz inverse problem of isostasy based on generating the isostatic gravity disturbances which equal zero.The formulation of the problem is done under the assumption of varying Moho depths T while adopting a constant value of the Moho density contrast 

Figure 1 .
Figure 1.The gravity disturbances computed globally on a 1 × 1 arc-deg surface grid using the EGM08 coefficients complete to degree 180 of spherical harmonics.Values are in mGal.

Figure 2 .
Figure 2. The topographic gravity correction computed globally on a 1 × 1 arc-deg surface grid using the DTM2006.0 coefficients complete to degree/order 180.Values are in mGal.

Figure 3 .
Figure 3.The bathymetric stripping gravity correction computed globally on a 1 × 1 arc-deg surface grid using the DTM2006.0 coefficients complete to degree/order 180.Values are in mGal.

Figure 4 .
Figure 4.The ice density contrast stripping gravity correction computed globally on a 1 × 1 arc-deg surface grid using the DTM2006.0 topographic and ice-thickness coefficients complete to degree/order 180.Values are in mGal.

Figure 5 .
Figure 5.The soft (a) and hard (b) sediments density contrast stripping gravity corrections computed globally on a 1 × 1 arc-deg surface grid using the coefficients of global sediment model complete to degree/order of 90 generated from the CRUST2.0density and thickness data of soft and hard sediments.Values are in mGal.

Figure 6 .
Figure 6.The refined gravity disturbances obtained after applying the topographic and stripping gravity corrections due to the ocean, ice and sediment density contrasts.The gravity data were computed on a 1 × 1 arc-deg grid of surface points.Values are in mGal.

Figure 7 .Figure 8 .
Figure 7.The Moho depths determined based on solving the VMM inverse problem of isostasy using the gravity disturbances corrected for the gravitational contributions of topography and density contrasts due to ocean, ice and sediments.Values are in km.

Figure 9 .
Figure 9.The Moho depth differences computed using the isostatic gravity disturbances and corresponding gravity anomalies.Values are in m.