The gravimetric contribution to the Moho estimation in the presence of vertical density variations

The Moho surface, namely the density discontinuity between crust and mantle, is traditionally studied by seismic methods. However, gravity information can contribute to the Moho estimation and, more generally, to the crustal modeling. The contribution is twofold. First, gravimetry generally provides observations with much lower errors than those implied by the mass density uncertainty and other geophysical assumptions. This means that it can be used to validate existing Moho and/ or crustal models by forward modeling. Second, gravity inversion is able to provide diffused (not localized) information on the mass distribution, both regionally and globally (thanks to dedicated satellite gravity missions). However, this information is weak due to its intrinsic ill-posedness. This means that it can be used to correct and spatially interpolate existing models, and to complement seismic, magnetic and geological information to create new models. In this work, the problem of estimating the Moho surface by gravity inversion assuming a two-layer model with lateral and vertical density variations is treated at a regional level. The approach consists in linearizing the forward modeling around a reference Moho at a constant depth and then inverting it through a Wiener filter. This is standard in case of two layers with homogeneous density distributions (or with lateral density variations), while it requires some additional considerations and algorithm modifications in case of vertical density variations. The basic idea is to “condensate” the masses inside the Moho undulation on the reference surface used for the linearization, thus requiring the setup of an iterative procedure. A strategy to introduce seismic information into this inversion procedure is proposed too, with the aim of improving the a priori density modeling. A closed loop test is presented for the algorithm assessment, showing the improvement with respect to a standard approach and the capability of the proposed algorithm to reconstruct the originally simulated Moho undulation by also fitting the gravity and seismic data at a level that is consistent with their observation noise.

. During this period, Gutenberg advanced his hypothesis about mountain building and its relationship with the deep structures of the Earth (see, e.g. Gutenberg and Richter 1954). At the same time, the integration between geology and geophysics was taking place because of the need in specific exploration problems, mainly related to the hydrocarbon industry, within the seismic prospecting community. After this period, the development of crustal studies was associated with the appearance and wide use of the so called "deep seismic sounding". This technique allowed the determination of averaged crustal structure over large regions, the construction of seismic sections along specifically located seismic profiles, and finally in 1982 the estimation of a first global Moho map (Soller et al. 1982). After that a noteworthy result, at a global scale, was the CRUST5.1 model (Mooney et al. 1998) based on seismic refraction and, later on, its updated versions, i.e. CRUST2.0 (Bassin et al. 2000) and CRUST1.0 (Laske et al. 2013). These models describe the crust structure by giving information on thickness and density of a number of global components (e.g. ice, oceans, soft and hard sediments, upper, middle, and lower crust) on a grid with a resolution of 5 • , 2 • , and 1 • , respectively. Moreover, for each cell of the grid, the crustal type (e.g. oceanic, continental plateaus, rifts, orogenetic regions, etc.) and the upper mantle density are also given. The models of the CRUST series are based on seismic refraction data published from 1950, on a detailed compilation of ice and sediment thickness, and on statistical predictions for regions such as most of Africa, South America, Greenland, and oceans where no or very few seismic measurements were available. The models are delivered without any error estimate, except the general assumption that seismic observation errors are approximately equal to the 10% of the Moho depth (Christensen and Mooney 1995). The problem of the Moho retrieval from seismic methods is well revisited, amongst many others, for instance in Carbonell et al. (2013). Parallel to seismic or seismologic investigations, volcanologists, petrologists, and mineralogists contributed to the study of the problem. From these studies, the crust was assumed to have compressional velocities V p of about 5.5-6.5 km/s and a density of about 2700 kg/m 3 , while the mantle presents V p -velocities of about 7.8 km/s or even more and a density of 3300-3400 kg/m 3 (Meissner 1973).
Also gravity observations have been exploited to get information on the depth and structure of the Moho. Again, apart from previous studies related to the problem of isostasy, which date back to 1850 with the theories of Herschel, Pratt and Airy and the subsequent refinements due to Hayford and Heiskanen (see Watts 2001 for details), the first attempts to determine the Moho depth, mainly from gravity anomalies, date back to 1950 thanks to the work of Heiskanen (1953). In 1960, Talwani derived rapid formulas to compute the gravitational effect due to bodies of arbitrary shape and applied them to estimate the crustal structure of the Earth by means of a trial and error method (Talwani et al. 1959;Talwani and Ewing 1960). A step forward have been done by Oldenburg (1974) who, rearranging the formula developed by Parker (1973) to compute the gravitational anomaly caused by a two-dimensional uneven layer of material, derived an iterative procedure based on Fast Fourier Transform to estimate the shape of the perturbing body given the anomaly. In the 90s, the Parker-Oldenburg forward and inverse modeling has been performed in three-dimensional (3D) space with large datasets. In the same period, several algorithms were developed with the aim of retrieving the mass density distribution from gravity anomalies, instead of the geometry of a surface separating two layers (e.g. Li and Oldenburg 1998). These methods are usually classified as linear inverse problems, due to the linear relation between the density and the gravity signal once the geometry is known (Blakely 1996), and rely on the discretization of the density distribution by elementary bodies, such as point masses, prisms, tesseroids, etc. They had very important improvements, going from the application of several regularizations (Boulanger and Chouteau 2001) to the more complex integration with a-priori geological information in a Bayesian framework (Marchetti et al. 2019;Reguzzoni et al. 2019;Rossi et al. 2015;Lane et al. 2007). However, here we will not focus on these linear methods since they are usually adopted in exploration geophysics to model complex crustal structures and not just to estimate the Moho geometry. Always in the 90s several other methods started to be proposed in the literature, for instance, Braitenberg et al. (1997) suggested an iterative method based on an isostatic assumptions, and Barzaghi et al. (1992) proposed an inversion method based on a Least Squares Collocation approach. Later, Strang van Hees (2000) derived a linearized expression relating the coefficients of a spherical harmonic expansion of the anomalous potential to those of the Moho depth. This allows for global inversions (see, for instance, Bagherbandi 2011, 2017;Reguzzoni and Sampietro 2015;Reguzzoni et al. 2013). Reguzzoni and Sampietro (2010) generalized the method proposed by Barzaghi et al. (1992) by extending it from the inversion of gravity anomalies to other functionals of the gravitational potential (such as the potential itself and its first and second vertical derivatives) and by applying it in the frequency domain.
Apart from the linear methods, which, as written before, aim at reconstructing the 3D mass density distribution, all the other algorithms to recover the Moho depth from gravity observations can be basically traced back to the socalled two-layer model, i.e. to the estimation of the interface between two layers with different densities. If the density contrast between the two layers is supposed known and constant, the uniqueness of the solution can be theoretically guaranteed at both global and local scale (Barzaghi and Sansò 1988). The uniqueness can be also ensured in case of vertical density variations, under suitable conditions (Sampietro and Sansò 2012;Barzaghi and Biagi 2014). This basically means that if we are able to reduce the observed gravitational field to the effect of just two layers with a given density contrast there is only one geometry defining the boundary between the two layers, given the observed signal. Note that, even when dealing with such a simple case, the inverse problem is still exponentially ill-posed (Sansò et al. 2018) and its solution should be properly regularized.
The general procedure to estimate the Moho from gravity data can be, therefore, summarized as follows: one has first to reduce the gravitational observations using the most accurate a priori geological information (possibly gravity independent), thus simplifying the initial problem in a twolayer problem. After that, the residual field can be inverted by means of (linearized) equations relating the observed functional with the unknown Moho depth. Possible instabilities, due to observation errors, are managed by means of a proper regularization operator.
In this paper, we follow this scheme placing the focus on how to deal with the presence of vertical density variations. In particular, in Sect. 2, we present a standard approach to solve an inverse gravimetric problem by linearizing Newton's equation around a reference constant Moho discontinuity surface. The solution is typically given by assuming that the mass density is homogeneous inside the crust and mantle layer, or at most by considering some lateral variations. We propose two alternative strategies to generalize the solution, the former is the simpler and consists in evaluating the density at the reference Moho, the latter by averaging the density between the reference and the actual Moho. Since the Moho is unknown, the latter solution implies an iterative scheme that is explained in Sect. 2. In Sect. 3, the problem is generalized assuming that the density profiles along the vertical direction are not perfectly known and an iterative scheme to estimate both the corrections to the density profiles and the Moho depth is presented. The three proposed algorithms are then compared in closed-loop simulations to assess their performances. These simulations are presented in Sect. 4. Finally, some final conclusions and remarks on possible further generalizations are provided in Sect. 5.

Moho gravimetric inversion with density variations
First of all, we have to distinguish between local and global inversion problems, assuming to work with anomalous quantities because the contribution of the normal gravity field has been previously removed from the data. The former are typically solved in planar approximation and the problem is formalized in Cartesian coordinates (x, y, z) . The latter can be solved in spherical approximation, typically using spherical coordinates, namely latitude , longitude , and radial distance r. In both cases, the Moho discontinuity surface D can be modeled as the sum of a reference constant surface D and an undulation D , i.e. z MOHO = D(x, y) =D + D (x, y) or r MOHO = D( , ) =D + D( , ) in planar or spherical approximations, respectively. This means that in planar approximation D is a plane, while in spherical approximation it is a sphere. Let us consider the classical scheme of the so-called twolayer model in planar approximation (see Fig. 1). Let us concentrate on a parallelepiped-shaped volume containing only crystalline crust above the Moho and upper mantle below the Moho (Fig. 1a). Let us assume that the observed gravitational signal has been previously reduced by the contribution of topography, bathymetry, sediments, and any other masses outside the considered volume, neglecting any border effect due to possible mismodeling. We will not discuss the problems related to this reduction, called stripping, in this paper, but we will focus on the Moho estimation assuming that the gravitational signal has been correctly stripped. We will call T the anomalous potential generated by this volume, and we will consider the vertical gravitational acceleration T z as observed signal, namely the derivative of the anomalous potential T with respect to the vertical z direction in planar approximation. The mass anomaly due to the Moho undulation ( Fig. 1b) is obtained by subtracting the reference mass distribution with a reference constant Moho from the actual one. If the crust and mantle densities are known, the result of this subtraction is a mass anomaly only between the reference Moho D and the actual one. Analogously, the vertical gravitational acceleration produced by the Moho anomaly only, that we will call T z , can be computed as the difference between T z and the sum of the gravitational signal produced by the crustal parallelepiped above the reference Moho and the mantle parallelepiped below the reference Moho, i.e. defining F the forward operator where C and M are the mass density distribution of the crust and mantle, respectively. D + or D − denote whether the same forward operator F is applied to the parallelepiped above or below the reference Moho D , respectively. For the sake of notation simplicity, we will now define the density contrast as if the densities have both lateral and vertical variations. Now we will first write the forward operator that transforms the Moho undulation into its corresponding gravitational signal T z , and then we will invert the relation (with a proper regularization) to find the solution of the inverse gravimetric problem. Starting from the usual simplified scenario where crust and mantle densities are constant between the actual and the reference Moho, namely the density contrast is constant too, the Newton equation describing the forward modeling of the anomalous potential in planar approximation can be written as and, therefore, the forward operator of the gravity disturbances is easily obtained by deriving Eq. 2 with respect to z, i.e. (1) where A is the domain area of the study region, G is the gravitational constant, x, y and z are the coordinates of the observation points and , and are the coordinates of the running point inside the masses. Equation 3 can be linearized with respect to D around the reference Moho D (namely D = 0 ). Using the Leibnitz rule, one can write If the gravitational signal is computed at a constant height from the reference Moho, h = z −D , Eq. 4 can be seen as a 2D convolution of the Moho undulation with the following kernel Sampietro 2010, 2012) This forward modeling can be, therefore, performed by a 2D Fourier transform, significantly reducing the computational burden, namely where F(⋅) and F(⋅) −1 are the 2D Fourier transform and its inverse, respectively.
Let us now generalize the definition of the crust and mantle densities. If they depend on the planar x and y coordinates only, the forward modeling of the Moho undulation is not changing very much. By defining the new variable the gravitational signal can be computed again by Fourier transform as The problem is slightly more complicated if the crust and mantle densities have vertical variations. We will see at least two ways of dealing with such a modification. The former consists in applying the linearizazion of Eq. 3 to the density contrast too. This means to evaluate the density contrast at the reference Moho depth, namely The latter consists in assuming that the density contrast between the actual Moho and the reference one is constant in the vertical direction, e.g. it is equal to the mean value inside the Moho undulation, namely In both cases, the forward modeling of the Moho undulation can be reconducted to Eq. 8. However, the two approaches are different if one wants to invert this forward operator to estimate the Moho from a given gravitational signal. Before commenting their differences, let us first face the problem of deriving an inversion operator from the forward one. Since Eq. 8 in its general form is a linear expression relating the product of Moho depth and density with the gravitational signal, its inverse is straightforwardly derived, i.e.
where T obs z is obtained from the observed vertical gravitational acceleration T obs z using Eq. 1. Note that the observations have to be further reduced at a constant height. The estimated Moho depth is then derived as However, the inverse operator of Eq. 11 is strongly unstable and can significantly amplify the noise of the observed gravitational signal. To overcome this problem, a regularization is required. According to the Wiener-Kolmogorov principle, i.e. the minimization of the mean square estimation error, the optimal regularization can be found by implementing a Wiener filter and, therefore, using the following inversion operator where K denotes the 2D Fourier transform of the convolution kernel and f 1 , f 2 are the two frequency components. S is the power spectrum of the variable , i.e. the Fourier transform of its covariance function assuming stationarity, and S is the observation noise power spectrum. The signal spectrum is usually defined starting from some a priori information on the Moho and density regularity, while the noise spectrum mainly depends on the instruments for acquiring the gravitational observations. Note that from such a stochastic modeling, one can also derive the spectrum of the Moho estimation error e = D −D and the corresponding covariance function by applying inverse Fourier transform. Now we can go back to the differences between the two ways of dealing with vertical density variations. The proposed inversion operator in Eq. 13 is suitable only for the first solution, i.e. fixing the density contrast equal to the value at the reference Moho depth (see Eq. 9). In the other solution, the density contrast depends on the Moho depth: in the forward modeling the Moho depth is known, but in the gravimetric inversion it is unknown, and, therefore, the density contrast is unknown too (see Eq. 10). For this reason, an iterative procedure is proposed where the density contrast is evaluated from the Moho depth estimated at the previous iteration, namely where is the iteration index, = 1, 2, … , L . The procedure is initialized by considering the density contrast at the reference Moho depth, i.e. assuming D(x, y) = 0 . The procedure is stopped when the variation of the estimated Moho depth between two consecutive iterations is negligible with respect to the estimation error. Actually, both solutions could be further iterated by redefining the reference Moho depth as equal to the mean of the final Moho estimate. Again, these "external" iterations have to be repeated until convergence.
A comment is due about the relation between the preliminary reduction procedure to isolate the Moho gravitational signal and the adopted inversion procedure. In both solutions, they have to be consistent. In the former, this consistency is satisfied because the reduction is performed by considering the actual densities of crust and mantle (see Eq. 1), and the inversion by considering the density contrast at the reference Moho (see Eq. 9). In the latter, there is no consistency and the reduction has to be corrected to consider that the density contrast inside the Moho undulation is not constant along the vertical direction, as it is instead assumed in the inversion procedure according to the approximation in Eq. 10. Since the Moho undulation changes iteration by iteration, also the reduction step has to be corrected at any iteration. Recalling Eq. 1, this means that y) at the iteration is given in Eq. 15. D + ( −1) or D − ( −1) denote whether the forward is computed for the Moho undulation above or below the reference Moho D , respectively. Note that the Moho undulation is the one computed at the previous iteration ( − 1).
It has also to be underlined that the uniqueness property of the two-layer model is preserved in the sense that, given the gravitational observations, there is just one and only one generating it, but of course there are many couples and D giving the same . Before concluding the section, we would like just to mention what happens in the global case in spherical approximation. The forward operator of the Moho undulation can now be expressed in terms of spherical harmonic coefficients, namely one can write the following linearized relation (Reguzzoni et al. 2013): where E is the mean Earth density, R E is its mean radius, T nm are the spherical harmonic coefficients of degree n and order m of the anomalous potential and nm are the spherical harmonic coefficient of the product between the Moho depth D and the density contrast. To derive this expression it is required that the density contrast does not depend on the radial direction between the reference Moho and the actual one. For this reason, if the density has a radial variation, the mean density contrast is considered in the Moho undulation, which is analogous to the second approximation proposed in the planar scenario, see Eq. 10. Also in the global case, an iterative solution is feasible, and this is actually the core of the solution proposed for the computation of the GEMMA Moho model .
Finally, we would like to recall the widespread Parker-Oldenburg inversion scheme, for the sake of comparison with the proposed linearized solutions. This approach was devised to deal with density distributions with at most lateral variations and, therefore, approximations like those in Eqs. 9 and 10 are required to manage vertical variations too. According to the Parker-Oldenburg modeling (Parker 1973), the forward operator of the Moho undulation for the gravity disturbances in planar approximation can be written as where f is the wave number, i.e. the modulus of the frequency vector f = f 1 , f 2 . The corresponding inversion operator (Oldenburg 1974) results and again this relation can be solved by means of an iterative scheme, updating the Moho undulation at each iteration (Gómez-Ortiz and Agarwal 2005).

Seismic data integration into the Moho gravimetric inversion
In this section, we would like to face the problem of integrating some seismic points in the Moho inversion schemes previously presented. Basically, a seismic point can be seen as a direct observation of the Moho depth. Since from the gravimetric inversion we are able to estimate the product between the Moho depth and the density contrast, the additional knowledge of the Moho depth in some points, in combination with gravimetric data, provides an indirect information on the density contrast. In the discussion of Sect. 2, we assumed that the density profiles along the vertical direction are perfectly known, while actually they are not. We could, therefore, use the additional seismic data to assess the given density profiles and possibly calibrate them according to some parameters, for example, biases and scale factors. In a simplified scenario, we can assume that the upper mantle close to the Moho surface has negligible vertical density variations, while the crystalline crust above the Moho surface can be divided into a set of "geological provinces", i.e. regions where all points have the same density profile along the vertical direction. Examples of such geological provinces and the corresponding functions relating the mass density with the depth can be found in Mooney et al. (1998). Since these functions are just approximate, they could be "calibrated" by introducing two parameters for each of them, namely a scale factor and a bias. According to these hypotheses, the crustal density can be modeled as follows: where i = 1, 2, … , N is the index identifying the geological provinces in the study area, i (x, y) is the characteristic function defining the domain of the geological province, C,i (z) is the given density profile and h i , k i are the corresponding scale factor and bias, respectively. Since any geological province is assumed to have a connected domain in a topological sense, different geological provinces of the same type can be present in the study area but, thanks to the calibration parameters, they can have slightly different density profiles along the vertical direction. In other words, it is possible that This crustal modeling has an obvious impact on the data reduction, that is slightly different for the two proposed inversion solutions, compare Eq. 1 with Eq. 16. Let us first consider the reduction in Eq. 1 and then generalize it to Eq. 16. Being the forward operator linear with respect to the density, the reduced gravitational acceleration due to the Moho anomaly only is given by where T M z is the reference upper mantle signal, T C,i z is the reference crust signal due to the ith geological province using its uncalibrated density profile, and T 1,i z is the same as T C,i z but using a unitary density profile, i.e. computed with C,i (z) = 1 . Now, since the inversion operator of Eq. 13 is linear with respect to the input signal, we can separately invert each term of the summation in Eq. 22, obtaining where ̂ , ̂ M , ̂ C,i and ̂ 1,i are computed from the inversion of T obs z , T M z , T C,i z and T 1,i z , respectively, using the operator of Eq. 13. Substituting Eqs. 21 and 23 into Eq. 12, we finally get the estimated Moho depth: This solution is suitable for the first strategy to manage vertical density variations, i.e. the one fixing the density contrast equal to the value at the reference Moho depth. As for the second strategy, i.e. the one considering the mean crustal . density between the reference Moho and the actual one, the data reduction has to be revised at each iteration , correcting Eqs. 1-16 through the additional term T ( ) z,corr . This is just a function of the scale factors h i for i = 1, 2, … , N , namely where By inverting T C,i,( ) z,corr , one can derive the corresponding ̂ C,i,( ) corr that has to be taken into account in the final inversion formula: Differently from Eq. 24, iterations are required to obtain the final solution from Eq. 27, because of the need of updating the data reduction through ̂ C,i,( ) corr and the mean crustal density ̄( ) C,i , both depending on the Moho undulation estimated at the previous iteration D ( −1) = � D ( −1) −D . The procedure can be initialized with D (0) = 0.
Since h i , k i are unknown for any i = 1, 2, … , N , the solution of Eq. 24 or Eq. 27 can be computed only if we suppose to know the Moho depth from external information, such as seismic data, at least in two points for each geological provinces in the study region. This solution can be obtained by least-squares adjustment, in case applying a Tikhonov regularization to keep the biases and scale factors as close as possible to 0 and 1, respectively. Once the parameters h i , k i have been estimated, they can be substituted into Eq. 24 or Eq. 27 to get the final interpolated Moho depth at all points (x, y) where the input gravity data were available. Note that if the forward operator of the Moho undulation is not linear with respect to , such as in the case of the Parker-Oldenburg algorithm, the parameters h i , k i have to be empirically estimated, e.g. by a trial-and-error technique. (25) .
Basically, the proposed algorithm allows to slightly change the crustal density profile in such a way that the estimated Moho is close to the seismic data and produces a gravimetric signal consistent with the observed one. In other words, this is a joint inversion of gravity and seismic data. However, since only few points are required for each geological province, we can state that it is a gravimetric inversion "weakly" combined with seismic data.

Numerical examples
The inversion schemes discussed in the previous sections were assessed by means of a closed-loop example. The test was performed by simulating the Moho surface over a 10 • × 10 • area in Central Europe, limited by 44.5 • N -54.5 • N and 5 • E -15 • E (see Fig. 2), on a regular grid with a horizontal resolution of 1∕8 • . The local Cartesian coordinates in planar approximation were defined as where , and h are the geodetic longitude, latitude and ellipsoidal height, respectively, ̄ and ̄ are the longitude and latitude of the center of the study area, and R is the radius of the local sphere.
The true geometry of the Moho depth was taken from the European Moho model (Grad et al. 2009), see Fig. 2a, the true crustal density was assumed to be linearly varying with depth according to three geological provinces ( i = 1, 2, 3 ), see Fig. 2b, c. Note that the first and the third geological provinces contain the same type of crust, meaning they have the same density profile, but they are independently processed because their domains are not connected. As for the mantle density, it was assumed to be homogeneously equal to 3300 kg/m 3 , without any lateral or vertical variation.
The aim of the closed-loop test is to show the performances of the proposed approaches by inverting the signal simulated from the model shown in Fig. 2 when providing different input data corresponding to a decreasing level of knowledge. In particular, the following three scenarios corresponding to different a priori density knowledge were tested: s1. density profile of each geological province equal to the true one; s2. density profile of each geological province corrupted by a scale factor h i = h (see Eq. 21), with the effect of underestimating the true density by 5%; s3. density profile of each geological province corrupted by a scale factor h i = h and a bias k i (see Eq. 21), with the effect of reducing the profiles slope by 5% and the density value at z = 0 by 2%.
The shape of these three a priori density profiles used as input to the inversion algorithm is shown in Fig. 3, while the corresponding h i and k i required to recover the true density profiles from the a-priori ones are reported in Table 1. The error introduced by the wrong a-priori density models, evaluated as the RMSE of the differences between the true and the a priori density profiles, is shown in Table 1 too. Moreover, we tested different versions of the algorithm to verify whether the increasing complexity of the gravity inversion procedure and the use of additional seismic data are really required. In particular, the following three versions were compared: . fixing an a priori density profile and a linearization depth D , and inverting the observed signal by means Eqs. 13 and 12; v2. fixing an a priori density profile and a linearization depth D , and inverting by iteratively applying Eqs. 15, 16, 17, 13 and 12 until convergence; v3. fixing a linearization depth D , and inverting by iteratively applying Eq. 27 to jointly estimate the Moho undulation ( D ) and the corrections ( h i , k i ) to the a-priori crustal density model, according to the input seismic information. The latter was simulated from the true D by adding a 1 km white noise.
All the three scenarios were tested with the three versions of the algorithm. It is worth to remark that we always inverted the same gravity signal, simulated at 1 km altitude from the true density profiles (see Fig. 2c) by adding a 5 mGal white noise as measurement error, and only changing the a priori knowledge on the density profiles. The grid of observations has the same extension and resolution of the Moho grid in order to apply the Discrete Fourier Transform in the inversion algorithms. The forward modeling to simulate the observed signal and reduce the data at each iteration of the algorithm (see Eqs. 15 and 16) was computed by summing the contributions of a set of rectangular prisms (Nagy et al. 2000) centered at the grid knots. Each prism has a size of 1∕8 • , corresponding to about 14 km×9 km in planar approximation. To introduce the vertical density variations, the profiles were discretized with a step of 0.1 km and a constant density was assigned to each prism in the vertical direction. The linearization depth D was chosen around the mean of the true Moho depth (see Fig. 2a) over the investigated region and fixed at D = −34 km.
Before commenting the results, a remark is due about the approximation introduced by the linearized forward operator (Eq. 4). In particular, we checked the impact of this approximation by comparing the anomalous signal T z (Eq. 16) due to the true Moho undulation D computed by the prism forward modeling and by the linearized formula. In the former case we applied the prism forward by considering the true density of crust and mantle, while in the latter case we applied Eq. 4 by introducing the mean density contrast ̃ (x, y) according to Eq. 10. Their differences, shown in Fig. 4, have a mean of 0.9 mGal , a standard deviation of  3.2 mGal and a minimum and a maximum of −4.3 mGal and 20.9 mGal , respectively. The errors are correlated with the Moho undulation, showing a worse behaviour of the approximated formula when increasing the undulation magnitude. This effect is due to the fact that, in practice, the linearized formula (Eqs. 4 and 10) corresponds to computing the effect due to the mass of the anomaly "condensed" at D . This difference has a systematic impact of the same magnitude of the noise of the simulated observations. However, from a practical point of view, this error is smaller than other sources of error that typically affect the gravimetric inversion, e.g. the data reduction for the effect of terrain, bathymetry, etc. .
To evaluate the results obtained with the different versions of the algorithm tested on the three possible scenarios, the mean, the standard deviation, the RMSE of the gravity residuals and the RMSE of the differences between the estimated and the true Moho depth were computed and reported in Table 2, while the estimated h i , k i and the RMSE of the differences between the estimated and true density profiles   with algorithm v3 are presented in Table 3. These results show that the algorithm is able to recover the Moho depth with an error of about 1 km . This error is fully compatible with the 5 mGal observation noise. However, this level of fitting of the true Moho depth can be reached only if additional seismic information is introduced to correct the a-priori density profiles through the estimation of the parameters h i and k i . Therefore, it is necessary to apply version v3, except when the a priori density profiles correspond to the true one (s1) for which, obviously, version v2 is enough. However, knowing perfectly the crustal density profile is possible only in an ideal world, thus the closed loop test shows the fundamental role of seismic information when the gravity inversion is applied to a real dataset.
The gravity fitting for the best cases is around 6.4 mGal and it is obtained by version v3 of the algorithm. This error is compatible with the 5 mGal simulated observation accuracy, where the slightly higher value could be attributed to the smoothing introduced by the power spectrum S of the field. In fact, the stationarity property of the field could not be locally verified, even if it is a good approximation over all the investigated area.
To better understand the improvement given by the different versions of the algorithm, we analyse the results obtained with scenario s3, the most realistic one. The residuals between the true and estimated Moho depth are shown in Fig. 5  (c) Province i = 3

Fig. 6
Estimated density profiles for each geological province using algorithm v3. The solutions corresponding to the three scenarios are represented by different line styles: solid (s1), dashed (s2) and dotted (s3), while the true profile is represented by black solid line stated, i.e. in a realistic case version v3 of the algorithm is strongly recommended. Note that it contributes to reduce not only the error standard deviation but also its mean, i.e. to correct possible systematic errors in the estimate. Furthermore, the solutions obtained for the different scenarios are very similar one to the other (see Table 2). In other words, by changing the input a-priori density knowledge, version v3 of the algorithm is able to correctly retrieve the same solution. Actually, the estimates of h i and k i are quite different from their true values for all scenarios, as can be seen by comparing Tables 1 and 3. This problem is mainly due to the correlation between the estimates of h i and k i , especially when they are applied to a linear density profile as in this example, see Eq. 21. However, from Fig. 6, representing the estimated density profiles with version v3, it can be seen that the impact of these "wrong" values of h i and k i on the final estimates of the density profiles is quite negligible.

Conclusions
This work aimed at studying the inverse gravimetric problem for a two-layer model with vertical density variations, like in the case of Moho depth estimation. The developed algorithm is based on the linearization of the Newton integral in planar approximation and its inversion by invoking the Wiener-Kolmogorov principle and by implementing a Wiener filter in the frequency domain. Three versions of the algorithm were proposed. The first one is based on the pure linearization of the Newton integral, which means to consider the density contrast at the reference Moho depth. The second one is based on the idea of condensing the masses inside the Moho undulation on the reference surface to have a more representative density contrast. The third one is like the second one, introducing into the model a bias and a scale factor of the density profiles to be jointly estimated with the Moho depth by also relying on information carried by external sources, like seismic-based observations. We performed a closed-loop test by developing the three versions of the algorithm and testing them on three different scenarios in which we changed the input information in terms of a-priori crustal density profiles. The results of the test show the capability to retrieve the Moho depth with an error of about 1 km compared to the true Moho shape used to simulate the signal, also when the a priori density profile is different from the true one, as usually happens when performing gravity inversion on real data. In this most general case, the closed loop test confirms that a joint inversion of gravimetric and seismic data is necessary when both the geometry and the density of the model have to be estimated.
In conclusion, the proposed algorithm for the Moho estimation to deal with vertical density variations seems to provide a numerical efficient solution to the problem, when working in the simplified scenario of a two layer model in planar approximation. The case of spherical approximation at a global level was already studied in the framework of the GEMMA project and it will be assessed in a closed-loop test like the one presented here in a future investigation.