The collocation approach to Moho estimate

In this paper, the collocation approach to Moho estimate is presented. This method is applied to the inversion of gravity data that can be complemented by Moho depth information coming from e.g. seismic information. In this context, a two layers model is considered and discussed in order to give a general theoretical framework for the inversion method. A body with two inner constant density layers and an inner separation surface between is considered and a uniqueness theorem is proved for the estimability of the separation surface given the gravity outside the body itself. Based on this result, a discussion is given on the estimation of the Moho depths based on terrestrial gravity observations. The observation equation is presented and its local planar approximation is derived. The application of the collocation method to the estimate of Moho depths is then studied and discussed in relationship to the planar observation equation. Also, numerical tests are presented. To this aim, the collocation inversion algorithm is implemented and tested on simulated data to prove its effectiveness. The results show that the proposed method is reliable provided that proper data reductions for model discrepancies are taken into account.


Introduction
The transition layer between the upper mantle and the lower part of the crust is defined as the Mohorovicic discontinuity (Moho) and ideally considered as a surface [Turcotte and Schubert 1982].On a global scale, this surface has sharp variations, ranging from over 60 km to less than 5 km depth.The density variation occurring across the Moho between the upper mantle and the lower crust is generally set at 0.5 g/cm 3 , based on the assumption that the upper mantle density is 3.4 g/cm 3 and the lower crust density is 2.9 g/cm 3 [Anderson 1989].The Moho transition layer causes seismic wave refraction and reflection and the related gravity signal as measured on the Earth surface can have a standard deviation ranging from 50 to 100 mGal.Moho depths can be estimated using both seismic and gravimetric inversion methods [Parker 1972, Lebedev et al. 2013].Seismic methods can give more accurate estimates that are not however homogeneously distributed over the Earth.On the contrary, grav-ity is quite homogeneously distributed over the entire Earth and only relatively small areas are un-surveyed [Pavlis and Rapp 1990].Furthermore, satellite dedicated gravity missions [Reigber et al. 1999, Albertella et al. 2002, Tapley et al. 2004] have made available a large amount of data that can be profitably used in combination with ground based gravity data [Shin et al. 2007].Recently, an approach based on the integration of gravity and seismic information has been proposed by Eshagh et al. [2011].They devised a method based on a stochastic combination of seismic and gravity Moho models.In this paper, the gravity estimation of the Moho based on the collocation method is presented [Krarup 1969, Moritz 1989, Barzaghi et al. 1992] which also allows integration of gravity and seismic information available in the investigated area.The basic idea of this method is to propagate the covariance structure of the Moho depth to the covariance of the observed gravity.Given gravity observations and Moho depth information (e.g. from seismics), the Moho collocation estimate can be obtained as a linear combination of these data.Also, collocation allows the computation of the variance of the estimation error that can be used as an indication of the estimated depth precision.This approach is here devised for a two layers model and a uniqueness theorem is proved as a theoretical background for the proper application of the method (see Appendix).In this context, the collocation inversion method for the two layers model is presented and the impact of model discrepancies on the Moho estimate is discussed.Computation problems are addressed too.The method effectiveness is then tested using simulated gravity data.are separated by a unique surface.Particularly, we consider a spherical two layers body.We also assume that this spherical body has a radius R = 6380 km and that the surface separating the external layer from the inner volume is at a depth ranging from 5 to 60 km.In a very simplified way, this can be assumed to be the reference model for the Earth and the Moho.Given this model, a new method is developed to obtain a Moho estimate on regional areas, e.g.areas having a 6° × 6° extension.Furthermore, the method and the equations are given in planar approximation.
Gravity observations are supposed to be known on the surface of the sphere representing the Earth.In addition, Moho depths can be given at scattered points.The density of the outer layer t C and the inner density t N are assumed as known and independent from the radial component.So, t C and t N depend only on the spherical coordinate v, v = ({,m).This is a special case in the theorem proved in the Appendix and thus it still holds that the inner separation surface is unique.The potential generated from this body in a point P on the external surface is expressed by , where N N and N C are the Newtonian operator applied to the inner volume and the outer layer respectively.If we define a normal potential is the Newtonian operator extended to the whole body), we can define the anomalous potential T as , where we have put t = t C − t M .This definition implicitly poses the density of the normal potential equal to the mantle density.Other choices can be given and are available in literature, for example by putting the 'normal' density equal to the Earth mean density.However, in our application the gravity anomaly is derived by T and a different choice of density just affects the mean dg that in any case will not be used.Thus T(P) is given by [Moritz 1990, Sjoberg 2009] P is a point on the external surface, R P , v P are the radial and angular coordinates of P: hypothesizing a spherical body, R P = .r PQ = ( ] PQ is the angle between OP and OQ, O being the center of the sphere).The term is the depth of the separation surface between the two layers.The equation for gravity dg = −^T/^R is then Since we can express dg in the following way with Q R and Q R−H running over the outer sphere surface and the inner separation surface.The planar approximation of these equations can be obtained in a reference system having the Z axis along the plumb line in the midpoint O of the investigated area and the (X,Y) plane tangent to the outer sphere in O (see Figure 1).We define n' as the inner surface normal vector and n as the unit vector along the radial direction.If we name dS as the area element of the inner separation surface, we have dSn' • n.The vector n can be further decomposed into e z + de x , being e z its component along the Z axis and de x the component onto the (X,Y) plane.It follows that n' If we assume to consider areas having maximum width of about 6°, de x , which is the tangent of the angle between n and e z , has a maximum value that is 1/20 of e z .This implies that, at the first order, we can approximate n' • n n' • e Z , that is dSn' • e Z .The second term in the equation giving dg can be thus simplified accordingly This integral can be further simplified by assuming 1 (this is true with an error which is, at maximum 1%).In this way, it follows that The same holds for the first integral.Hence that is where .This is a non-linear integral equation between H(x,y) and dg.In order to get it in a form that is suitable for applying collocation, we have to linearize it.To this aim, we define (E[+] is the mean operator) and we set H(x,y) = .Moreover, let assume that : note that in some areas this is not true and the resulting approximation could represent one possible limitation of this approach that must be carefully considered in its application.In any case it is needed for getting dg(x,y) as a linear functional of f(x,y) as required for applying collocation: more investigations on the approximation consequences will be discussed in the tests of Section 4. Therefore, considering only the linear terms in f, we obtain If we further consider the density contrast between the crust and the upper mantle, we set and .Thus, the equa-tion for dg(x,y) can be written as with In its linearized form, the dg(x,y) observed gravity is given as the sum of four terms that will be analyzed in the following.Integrals I 2 , I 3 are linear terms in f and dt while I 4 contains the product f(p,i)dt (p,i).This last integral is thus a second order term which is likely to be smaller than I 2 .So, in a first approximation, we assume that I 2 + I 4 I 2 .If we now disregard I 4 , we have In our hypotheses on f and dt, it holds that , which implies that .So, in principle, one can estimate from observed gravity.This is indeed possible only in theory but not in practice due to the standard definition of observed gravity anomalies (there is a mismatch in the reference density between our definition and the standard one based on the normal ellipsoidal field).Thus we assume to have an estimate of as given independently (e.g. from seismic information).
The last term to be discussed is I 3 .It gives the gravity effect due to density variations over the thickness , ,  ( , ), [ ] 0 layer.The impact of this term can be large as it can be assumed that, in the mean over the whole Earth, 30 km.Hence, this term must be locally modeled based on available crustal density variation data in the investigated area.
The amplitudes of the four integrals terms and the impact on the estimated f values will be further discussed in the numerical simulations given in paragraph five.

The collocation estimate method
It is well known [Menke 1989, Tarantola 2005] that the gravity inversion is unstable.Thus the estimation algorithm to be applied should be stable and smooth with respect to high frequency gravity signal components.Collocation method can be successfully applied to this aim.It is a stochastic filtering/predicting method that filters out high frequencies based on the covariance structure of the data.It has been widely applied in Geodesy [Tscherning 1985, Moritz 1989] particularly in view of the use of heterogeneous data.In our application, this implies that dg data can be used together with f known values (coming from e.g.seismic) to get an integrated estimate of f.This method can be properly applied to the observation equation .Thus, in order to get unbiased estimates, we need to model and remove the signals related to the I 3 and the I 4 terms that account for density variations.As it will be proved in the simulations, these terms give a significant signal having also a frequency signature that couples with I 2 .However, it can be reasonably assumed that I 4 ≪ I 2 and I 4 ≪ I 3 .Thus, the proposed estimation method is performed according to the following steps.
1) Estimates of and dt(x,y) are given using available information on t(x,y).
2) On observation points, the I 3 term effect is estimated based on the given values of dt(x,y) and (this last term can be derived from known geological information and/or from seismic profiles).This effect can be straightforwardly computed using the closed formulas of the gravity effect of a prism [MacMillan 1958].
3) Reduced gravity data are estimated as and, assuming that dg I (x,y,0) are estimated by collocation.
4) Using the estimated values, I 4 is computed and subtracted from dg I (x,y,0) to get dg II (x,y,0) = dg I (x,y,0) − .5) Then, the last computation step is accomplished using dg II (x,y,0) as input data and getting, by collocation, . Thus, in the end, In order to see how collocation can be applied to get the estimate of f, a brief presentation is given in the following.A more detailed description can be found in Barzaghi et al. [1992], Biagi [1998] and Forsberg [1984].
We assume that the following linear relationship between gravity (dg) and depth (f) holds Assuming that and are known, the collocation estimate of f can be obtained based on dg observations and f values given in some scattered points.
To get these estimates, we further have to assume that: 1) f is a weak stationary stochastic process, ergodic in the mean and in the covariance 2) the noises in gravity and depth, n g and n f are uncorrelated zero mean signals 3) the cross-correlation between signals and noises is zero.
Auto and cross covariances can then be summarized as follows Using these covariances and observations plus noise we can get the following auto-covariance matrices of the observations The collocation estimate of f [Moritz 1980, Barzaghi et al. 1992] in a point P k = (x k ,y k ) is given by with (3.9) To compute this estimate, we first have to estimate the empirical covariance function of dg which must be then modeled with proper positive definite model functions [Moritz 1980].Also, models for auto and cross-covariances must be defined in order to compute the matrices defining the f estimator.
This can be done via covariance propagation law applied to linear functionals [Moritz 1980].In our approach we assumed where L P [f] indicates the linear functional relating f to dg.

It follows by covariance propagation law that
A suitable model for C ff can be (J 0 : zero order Bessel function) which implies [Biagi 1998] Another model that can be used for the auto-covariance of dg is (J 1 : first order Bessel function) which gives The A and the a values contained in these models are tuned to fit the empirical estimated covariance values of dg.As stated above, models for C ff and C fdg can be then derived, allowing the computation of the matrices contained in the estimation formula for f.

Some comments on the efficient computation of the estimation formula
One numerical major problem in computing the collocation estimate of f is the inversion of the covariance matrix of the data C ll .Some standard data configurations are discussed here.We firstly assume that the number of known Moho depths N H coming from seismic is at maximum 200.We then consider three different gravity data distributions: 1.A quite poor gravity data coverage and depth information: 2. High gravity data coverage and no depth information: 3. High gravity data coverage and depth information: Case 1. Can be handled using standard techniques such as Cholevsky decomposition [see e.g.Press et al. 1992].The efficient numerical computations of case 2 and 3 can be obtained if we assume that gravity data are regularly gridded.In this case, the covariance matrix has a Toeplitz/Toeplitz structure: based on this particular structure, in case 2, Fast Collocation can be applied [Bottoni and Barzaghi 1993].Case 3 must be discussed in more details.Also in this case, we assume to have gridded gravity data.On the other hand, depth information are assumed to be given on scattered points.Under these assumptions, an efficient solution can be obtained by partitioning the solution.We have to solve the following system The partitioned solution is given by with In this way, the term m can be computed via Fast Collocation being dg on a grid.The same holds in com-, , , , where c i is the i-th column of , i = 1,..., N H .The memory allocation for this method is equal to (N g + N H )(1 + 2N H ) since we have to store the first row of the Toeplitz/Toeplitz matrix C dgdg , the full C ff and C fdg matrices, the C and X matrices and the vector m.This memory occupation is smaller than those implied by the full C ll matrix storage.Furthermore, in terms of computation effort we have (N H + 1) Fast Collocation solutions having N g dimension and 2 Cholesvky based solutions having N H dimension.The computation time is thus proportional to (N H + 1)N d 2 g logN dg + 2N H 3 which is less than the one require for solving the full C ll system by standard Cholevsky method.
A last remark is in order.Usually collocation is a stepwise procedure since it is an adaptive filter which is tuned by the empirical covariance.Thus, in the first step, the low frequency component is dominant while in the further steps the high frequencies become dominant.At each step, residuals on observation points are obtained as the difference between observed and collocation predicted values.The iterative process stops when residuals (e.g.gravity residuals) are uncorrelated.

The simulated tests
The devised method has been tested via numerical simulations that aim at verifying: -the properties and the stability of the inversion procedure under model/data consistency, -the impact of known H OBS data in the inversion, i.e. the improvements that one can get with respect to the inversion based on gravity data only, -the biases induced in the depth estimates when model inconsistencies are present, like for example, errors in the gravity contrast model or in the mean Moho depth.
A surface depth has been simulated on a regular (x,y) grid with 10 km grid knot over a rectangular area 600 km × 400 km.The depth values are given with respect to mean depth while the f values have been generated according to the procedure described in Barzaghi et al. [1992].The a-priori covariance structure of the f values has been fixed according to the following parameters The simulated depths (Figure 2) have the following statistics v(f) = 6.4 km min(H)= −43.67 km max(H)=−16.00km As one can see in Figure 2, the resulting depth surface is smooth and regular which seems to be quite far from the real Moho behavior, at least in some particular areas where high frequencies features are present (e.g. the Alpine region).However, this is not so relevant since in the devised tests we are aiming at defining the impact of model inconsistencies on the inversion procedure.This can be critically evaluated even if high frequency pattern is not present in the target surface.Generally, the simulated depth does not respect the condition adopted for the (2.9) linearization, i.e.
: this will allow a sensitivity assessment of the proposed al- Starting from the simulated depths, the gravimetric effect has been computed according to two different density contrast models.In the first case a constant density contrast t(x,y) = t 0 = −0.52 g/cm 3 has been considered.In the second case, it is assumed that the density contrast varies linearly with x: In both the cases (Figure 3), the gravimetric effect has been computed at z = 0 by prism integration for-mula [MacMillan 1958] over the same grid used for the depth values.In the following, the two computed gravity fields will be named G C (x,y) (constant t 0 value) and G L (x,y) (linearly varying t L value).Furthermore, three depth sections have been computed based on the H(x,y) grid.In this way, we simulated 2501 H(x,y) grid values, 2501 gravity grid values and a total of 100 scattered depth observations on the three distinct sections.
These preliminary simulations are aimed at the discussion of the theoretical model and the analysis of the implemented software: therefore, no noise has been added to the simulated observations.Three different inversions have been carried out: 1) Integrated inversion based on G C (x,y) gravity , ; .g/cm , .g/cm ; .g/cm and known depth values assuming a constant gravity model t 0 .This is the case when data and density are consistent.
2) Inversion based on G C (x,y) gravity data only assuming the constant gravity model t 0 .Also in this case data and density are consistent.
3) Integrated inversion based on G L (x,y) gravity and known depth values assuming that density is linearly varying following the t L (x,y) rule.Once again data and gravity model are consistent.
4) Integrated inversion based on G L (x,y) gravity and known depth values assuming that density is constant and equal to t 0 .This is the critical case where we have a mismatch between data and density model.
In all the inversions, the model covariance function of the gravity field has been fitted on the empirical covariance function of the simulated gravity observations [Barzaghi and Sansò 1983]: therefore, it has been propagated by (3.14) and (3.16) to model the cross covariance between gravity and depth and the covariance of depth.Data inversions 1) and 2) are made for testing the impact of known depths in the collocation inversion procedure.Test 3) is performed to check for the effectiveness of the iterative inversion procedure described in paragraph 4, assuming to have a known density model consistent with the available data.Finally, test 4) is set to estimate the impact of improper density information on the integrated estimated depths.
Inversion 1) should provide the benchmark results, for three reasons.The adopted density model in the in-version is consistent with the simulated model, all the simulated observations are used and the constant density model does not require either approximations or iterations to manage I 3 and I 4 terms.
In all the inversion procedures, as it would be in a real case, the mean depth has been estimated using the 100 given depth values.The obtained value is .In this way, the inversions are based on a biased mean depth.So, the impact of this bias on the results can be investigated too.
Furthermore, as a general remark on all the inversion procedures, it must be underlined that two iteration steps were necessary to get the final estimates since first step residuals were still spatially correlated (which was not the case for the second step residuals).
The estimated depths have been then compared with the simulated ones in order to evaluate their precisions and accuracies (to avoid edge effects that can give biased statistics, comparisons are carried out in the 500 × 300 inner area).Statistics of the differences for all the four inversion procedures are listed in Table 1.Statistics over the 100 known depths are also listed when such data are used in the different inversion procedures.Figures 4 to 6 report the results provided by inversion 3, that are the most interesting.
The following remarks are in order.The results obtained for inversion 1 (see row one in the a-priori bias remains in the final estimates, because no depth observations are given to tune the mean.In this case the mean of the differences is 2.2 km (see row three   in Table 1).However, it must be stressed that the f values are correctly estimated as it results from the standard deviation which is equal to the one obtained in inversion 1.The iterative procedure used in inversion 3 to correct for density heterogeneities leads to satisfactory results.It must be remarked that the computation and correction for the I 4 term improves the results, even if this term is remarkably smaller than I 3 (Table 1, row four are statistics base on modeling I 3 term only; in row five, I 3 and I 4 are used in the estimation procedure).
Generally, in cases where no model error is present in the inversion (inversions from 1 to 3) the standard deviation of the residuals (0.4-0.5 km) is less than 10 % of the standard deviation of the simulated depths (6.4 km); moreover, the approximation adopted to linearize the observation equations does not affect significantly the results: also on this regard, results are completely satisfying.
The adoption of an improper model for density heterogeneities has a significant impact on the estimated depths (see Table 1, row seven).This is true even in the case of the integrated inversion based on gravity and depth observations.The impact of density model mismatch can be evaluated in the residuals computed on the observed depths that show high discrepancies.Hence, the observed depth values cannot compensate for an improper density model thus confirming the importance of a reliable correction for existing density heterogeneities (see also row eight in Table 1; even the residuals on given depth points are worse than in the previous cases).

Conclusions
The collocation method can be profitably used to combine gravity and depth information in order to get reliable Moho estimate.In this framework, a theoretical background has been provided through a uniqueness theorem for a two layers body.Under quite general hypotheses on the body inner density distribution, it can be proved that the separation surface between the two layers can be uniquely estimated.Based on this theorem, the estimation of the Moho is discussed and the equations relating gravity and Moho depth are derived in planar approximation.In this context, a procedure to Moho estimate based on collocation approach is presented.This iterative approach has been then tested via numerical simulation.A simulated Moho is computed and two density models are assumed to compute the direct gravity effect to be used in the inversion procedures.These tests gave relevant insight in the devised inversion method.Reliable results have been obtained when density models are consistent with gravity simulated data.It can be thus assumed that, under unbiased conditions, the method is stable and can be applied to Moho depth estimate.Also, the impact of depth information proved to be relevant.Collocation allows easily the integration of different data types and this can improve the solution.Particularly, this applies to the bias in the mean depth that is reduced if depth data are considered.On the contrary, the solution based on gravity only cannot compensate for this bias.Another important remark is related to density information.It was clearly proved that density plays an important role in the inversion.If proper density information are available, the iterative inversion procedure based on collocation can give reliable results.On the other hand, when assuming inconsistency between density and simulated gravity data, we obtained remarkable distortions in the estimated depths.Furthermore, in this case, depth information are not effective in reducing the bias and even on the given depth data a significant discrepancy is present.
Thus, it can be concluded that the collocation based procedure can be applied for getting Moho depth estimates from gravity and depth information (from e.g.seismic).However, as it is for any other inversion method, improper density information can induce serious biases in the estimates even if depth data are included in the inversion.Two possible evolution lines of the collocation inversion method can be devised.One line is to investigate its behavior in local areas where high frequency Moho patterns are present.Being collocation quite a stable and regularizing technique [Barzaghi et al. 1992], this doesn't seem to be a too critical issue.The second possible application of this methodology is on global Moho estimates based on gravity data coming from the gravity dedicated satellite missions.Since proper spherical model equations can be quite easily derived and linearized [see e.g.Eshagh et al. 2011], global applications of the collocation method seem to be feasible.

Appendix. A uniqueness theorem for the two layers model
The gravity inversion problem can be defined as follows: to estimate the inner density distribution of a body given gravity observations in the spatial domain outside the body.As it is well known, this problem is an ill-posed problem [Moritz 1990] because uniqueness is not, in general, guaranteed.As an example, concentric spheres having the same total mass homogeneously distributed can vary their densities and volumes to give the same potential in the space outside the larger considered sphere.Thus, by inverting the gravity potential no unique solution can be estimated.To have uniqueness in the gravity inversion problem one has to impose proper regularizations and/or restrictions to the required solution [see e.g.Sansò 1980, Moritz 1990, Sjöberg 2009].In the framework of the gravimetric Moho estimation problem, a uniqueness theorem for the two layers model is given in the following.This theorem can be considered as an extension to a more general case of the one proved in Barzaghi and Sansò [1988].It can be used as a reasonable model for the gravity inversion to Moho estimate in the real case.
Two bodies B 1 and B 2 are given with two layers inner mass distributions and the same bounding surface S. The two inner volumes are B i1 and B i2 and the two crust layers are B e1 and B e2 .The two separation surfaces between the inner volumes and the crust layers S i1 and S i2 are considered as regular and are parameterized as S ik = R ik (v), with k = 1,2 and v = ({,m) (see Figure A.1).Let t i1 and t i2 be the inner part densities and t e1 and t e2 the crust densities of the two bodies.These densities are given as a function of (r,v).Assume that two functions t i (r,v), t e (r,v) are defined inside S and satisfy the following conditions Moreover, the following conditions hold:  .3)this body is equal to zero since where N(t) is the Newtonian operator.
The density function t has the following properties: Furthermore, having by hypothesis that , it follows We now define the two surface S i and S e are the bounding surface of B. We can further define S e+ and S i+ as those parts of S e and S i that bound the domain in B where t is positive.Similarly, we define S e− and S i− as those parts of S e and S i that bound the domain in B where t is negative.Furthermore, we define v+ and v− as the projection of S e+ and S e− onto the unit sphere v.
The potential generated by B is harmonic down to S e since the mass of this body is inside this surface.As stated before, the potential of B is equal to zero outside S. Due to the maximum and minimum property of the harmonic function, it is also equal to zero outside S e .It follows that t is in the orthogonal complement of L H

B
We now define a function k(v) onto S e , satisfying the following condition where H v 1/2 is the Sobolev space ½ of the function having domain on the unit sphere.
Based on the hypothesis on R i1 and R i2 , we can also assume that R e and R i are regular surfaces.It follows that .Furthermore, in the hypothesis stated for k(v), there is a unique harmonic function h(v) defined in B that on S e has the values assigned to k(v).If , it follows that .Hence Being the following decomposition for I holds where in in Since, obviously, I 5 can be written in the following way It holds that I 5 = 0. We suppose now that R e (v) differs from R i (v).Under this hypothesis, we want to define the values assumed by the remaining integrals.Since h(p) is harmonic in B, it can assume its minimum or maximum values only on the boundary, i.e. onto S e .In the limit for l + (v) → 1, we obtain 0 ≤ l(v) ≤ 1, that is .It follows that By making use of the properties of ^r t − , we have that I 4 < 0 which implies I 3 − I 4 > 0. The I 2 term is positive.If we define I II as we have I 2 < I II which implies I 1 − I 2 > I 1 − I II .I II can be further be decomposed in the following way Hence, we obtain that where M + is the positive mass included in the domain where t is positive.It follows that If we now consider the limiting value for l, we can write the following Since we finally have I 1 − I II > 0. So, if we assume that R e (v) doesn't coincide with R i (v), it holds that I = I 1 − I 2 I − 5 > 0. However, I must be equal to zero.This implies that R e (v) and R i (v) must coincide, i.e. that the interface surfaces in the two bodies are the same.

Figure 3 .
Figure 3. Gravity field with constant density model (top) and differences of the gravity field with the linear density model (bottom).
THE COLLOCATION APPROACH TO MOHO ESTIMATE

Figure 4 .
Figure 4. Integrated inversion 3 (linear density model and consistent correction of I 3 and I 4 terms): prediction errors.

Table 1 )
proved that depths can be correctly estimated even if the a-priori mean depth is biased.If scattered depth data are available, they can correct the initial bias and the final mean

Table 1 .
Statistics of the residualsD = H simulated − .H estimated t discrepancy between estimated and simulated depths is 0.5 km, which is negligible if compared to the mean depth of 30 km.If gravity data only are used (inversion 2) Tapley, B.D., S. Bettapur, M. Watkins and C. Reigber  (2004).The gravity recovery and climate experiment: Mission overview and early results, Geophys.
© 2014 by the Istituto Nazionale di Geofisica e Vulcanologia.All rights reserved.
can be proved that S i1 and S i2 are the same surface, i.e. thatProofWe consider a body B having a density distribution which is obtained by subtracting the B 2 from the B 1 density in every point.Outside S, the potential implied by