An Empirical Method for the Optimal Setting of the Potential Fields Inverse Problem

The use of potential field methods for geophysical exploration purposes is nowadays quite common: these techniques consent to retrieve geological knowledge over extended regions and can give complementary information where other invasive or expensive techniques, such as seismic acquisitions, fail (e.g., in the recovery of geometries of geological horizons beneath a thick salt layer). Recent dedicated satellite gravity and magnetic missions, such as GRACE, GOCE and SWARM together with the exploitation of offshore satellite altimetry and airborne/shipborne surveys, have paved the way to the realization of a variety of global models, characterized by spatial resolutions of about 4 km (both for gravity anomaly and lithosphere magnetic anomalies) and high accuracy (about 3-5$3\text{-}5$ mGal and 20 nT). These models are a valuable source of information to study the geological evolution and characterization of the lithosphere structure, especially at a regional scale. In the present work, some preliminary technical aspects related to the use of these models to perform three‐dimensional inversion are discussed, thus defining an empirical but rigorous procedure to set up gravity and magnetic inversion. In particular, we address the questions whether the classical planar approximation is acceptable for regional inversions or if a spherical one is required. We also provide guidance for choosing the best gravity functional (e.g., gravity anomalies or second radial derivative of the anomalous potential) and the optimal sizing of the three‐dimensional volume area to be modelled depending on the specific target investigated. The application of the proposed methods to the Mediterranean case study is also presented.


INTRODUCTION
Potential field methods are powerful tools to recover fundamental information on the Earth's crust structure. Gravity and magnetic data have been collected worldwide with dedicated satellite missions, and, thanks to these data, large-scale crustal structures are nowadays quite well known. To fully exploit these global datasets, the development of physically based modelling techniques is required. Inferring the Earth's crust structure, depicting the boundary of the geological units at depth and the stratification of the crust, requires in fact integrated approaches reconciling all the available measurements not only at satellite altitude but also (when possible) at a near-surface regional scale. Inverse gravimetric/magnetic problems aim at reconstructing subsurface properties given the observed field above the Earth's surface.
Geophysical literature on the solution of inverse problems related to gravity and magnetic data is extremely vast; in our work, we will cite just few seminal publications and suggest the papers of (Nabighian et al., 2005) for a more in-depth review. In general, gravity and magnetic methods comprise algorithms to estimate the three-dimensional distribution of density or susceptibility (e.g., Y. Li andOldenburg 1998a, 1996, for the gravity and magnetic field, respectively), solutions aimed at retrieving the geometry of the surface separating two layers of a given known property (see, e.g., Oldenburg, 1974, then evolved in more general surface geometry inversion in Galley et al., 2020, and the reference therein) or to more complex approaches, formalized, for example, by means of Bayes theorem (Bosch et al., 2006). All the above inversions methodologies, before tackling the gravimetric/magnetic problem, require the definition of a mathematical model of the sources usually based on approximations and simplifications. In the current work, we focus on this preliminary phase delineating a strategy to • numerically evaluate whether the classical planar approximation is acceptable for regional inversions or if a spherical one is required, • give the guidance for choosing the best potential field method/gravity functional (e.g., gravity anomalies or second radial derivative of the anomalous potential) to study the given target and • optimally sizing the inverse problem.
The setting of the above points is usually left to the experience of the operator performing the inversion (like in, e.g., Y. Li and Oldenburg 1998b) or, even worse, to the available inversion software; however, it can have important consequences on the inversion results. For instance, the proper setting of the spatial resolution of the model of the sources can reduce the computational time and power required to perform an inversion of orders of magnitude, while a wrong setting of the mesh padding zone can introduce important perturbations in the inversion solution. The literature on this topic is almost absent, apart from studies aimed at estimating the difference between planar and spherical approximation for a specific scenario (e.g. Kuhn et al., 2009), to compare gravity and gravity gradients (X. Li, 2001) and to study the spatial resolution of a given inverse problem (see An, 2012, and the references therein). To be more specific, in Kuhn et al. (2009), the Bouguer anomalies computed with spherical and planar approximation over Australia are compared, finding an almost constant bias over areas with moderate elevation changes. In the second work, X. Li (2001) provides some (empirical) technical procedures to compare the vertical resolutions of a gravimeter and a gravity gradiometer, while An (2012) proposed an algorithm to determine the spatial resolution of a general inverse problem.
The current work follows a principle similar to the one used in X. Li (2001), based on a set of forward problems, but extending the complexity of the targets (from simple shapes like spheres or cylinders to more realistic geological horizons) and the number of investigated questions. To fix the ideas and to simplify the discussion, we will focus on a specific inversion method, namely the Bayesian inversion presented in Sansó and Sampietro (2022), in which the investigated volume is discretized in a set of volumetric elements and the inversion is applied to recover the property of each element. However, the proposed approach is more general and can be applied to all the inversion methods identified in this section. Moreover, since the answers to the questions investigated in our work are not straightforward and are strongly related to the target (as emerged also from X. Li, 2001), the proposed approach is here applied to the Mediterranean Sea region with the dual objective to show a real numerical example and to define a benchmark for further similar studies.

METHOD
We present here briefly the inverse problem considered in the current work. We focus on problems at basin and regional scales (i.e., from about 50 km × 50 km to 3000 km × 3000 km or even more) trying to develop scalable solutions where the scope is the determination of the three-dimensional (3D) distribution of density and magnetic susceptibility inside a given volume. In the least-squares sense, the data misfit is computed by the classical quadratic form (see, e.g., Y. Li and Oldenburg, 1996): where is the vector of the observations, is the vector of unknowns, is the observation error covariance matrix and F is the general forward operator relating the physics of the sources to the available observations. In the case of gravity anomalies in planar approximation, the forward operator is given by (Blakely, 1996) where , and are the coordinates of the -element of the observations vector, is Newton's gravitational constant and , and are the integration variables running in the volume of the sources. For the magnetic field, supposing induced magnetization only, the forward operator is (Blakely, 1996) where is a constant depending on the used unit system, is the main inducing field vector and is the magnetic susceptibility. In principle, (2) and (3) should be closed integrals around all sources. However, when dealing with local and regional inversions, the considered volume should be limited both in the lateral extent and in depth. Moreover, ( , , ) and ( , , ) are usually discretized in a set of values and , each one representative of the mean density and susceptibility distributions inside an elementary volume. One of the most common discretization in planar approximation consists of dividing the whole volume in a set of rectangular prisms with size Δ × Δ × Δ .
Considering as observations in Equation (1) the gravity and magnetic fields, and assuming measurement errors for the two quantities independent with covariances and BB , respectively, Equation (1) becomes wherēis an × matrix, with number of gravity observations and number of discretized densities, containing at the element , the gravitational effect of the prism to the observation point (obtained by integrating ). Similarly, is an × matrix, obtained by integrating . Equation (4) represents the base of every joint inverse problem. At this point, we are left with a set of given observations , to be inverted in order to estimate a finite set of densities and magnetic susceptibilities , within the discretized study region . In such a problem, we can consider as known the observations as well as their accuracy and spatial resolution. On the contrary, the volume (defined by the coordinates of its vertices) and the size of the prisms in which it is discretized (namely Δ , Δ and Δ ) should be defined before the inversion. In addition, to avoid edge effects, an external border region of side Δ should also be considered.
An important remark has to be made on the forward operators defined in Equation (4). When dealing with such regional inversion problems, the gravitational effect of the Earth's normal field is removed (Moritz, 1990). When moving from the ellipsoidal to the planar approximation, the removal of the normal field can be approximated to the removal of an averaged constant field generated by an unknown 3D density model made by a set of Bouguer slabs (Sampietro, 2015).
At this point, two possibilities are left, either to work with density anomalies with respect to the unknown background model or to fit the gravitational field apart for a constant mean value (Sampietro & Capponi, 2019). Choosing the latter possibility, the first term of the data-misfit equation can be written as where E[ ] is the expectation operator applied to , is the number of observations, is the identity matrix of size and is an × all-one matrix. Basically, ( − )r epresents a zero-mean forward operator (Sansó & Sampietro, 2022). A similar argument can also be used for the magnetic field (at least for regions smaller than 2500 km × 2500 km). In fact, when trying to recover information on the Earth's crust from the magnetic field, the lithospheric contribution should be separated, among others, from the core one in the observed field. This is usually done by analysing the power spectrum (degree variances) of global magnetic field models expressed in terms of spherical harmonic series (Thébault et al., , 2021. Total magnetic internal field models show a break in the slope of degree variances between degrees 13 and 16 (Langel & Estes, 1982) roughly corresponding to wavelengths between 3000 and 2500 km (Backus et al., 1996). Longer wavelengths are interpreted as the core field and its secular variation, while shorter wavelengths are attributed to the lithospheric magnetic field. From a purely mathematical point of view, the difference with respect to what happens for the gravity field is that in the gravity case the effect of spherical harmonic coefficients for order m = 0 (rotational symmetry) and first even degree (2, 4, 6) are easily removed from the observations; while in the magnetic case, crustal field contributions from the coefficients up to degree 16 are not available since they are masked by the core field.
The size of the investigated volume, the regional extent of the measurements and their spectral content imply implicit assumptions on the sources that must be evaluated in order to set up the parameterization of the inversion problem. In our work, analysing the correct parameterization of the inverse problem, we will be led by a simple but rigorous principle: we will compare the 'observation error' and the 'model error' with the objective to define the initial model of the sources with a degree of complexity compatible with the observations' accuracy. In other words, the main idea is that if a certain simplification of the model entails a model error (in terms of standard deviation) smaller than the observation error, this means that the simplification can be safely applied because, in any case, its effect will be masked by the noise that affects the observations.
To do this comparison between observation and model errors, we consider the observation error as a fixed limit in the sense that, given the observations, the related error exists and cannot be modified or removed. So, to describe the proposed approach, it seems natural to, first of all, assess the accuracy of the potential field data and, after that, to evaluate the different model errors related to different possible simplifications and assumptions. This latter point is simulated by exploiting a realistic synthetic model of the considered volume and solving a set of forward problems by means of the ( − )ā nd ( − )̄operators. To evaluate these model errors, the effect of the studied volume (in terms of the grav/mag potential field) will be first computed by means of the forward operator in the most accurate way and then will be compared with the effect calculated by introducing one or more simplifications. Note that, once the gravity and magnetic effects of the synthetic model are computed with the above operators, the resultant fields will be zero mean thus the standard deviation calculation equals the root mean square one. We are therefore allowed to consider indifferently one of the two indexes while analysing the results. In the following, we will focus, for illustration purposes, on the Mediterranean Sea area as a test case.

THE MEDITERRANEAN SEA CASE STUDY
As a case study, we consider a region ranging between 6.2 • W to 36.3 • E and 29.7 • N to 46.1 • N covering the whole Mediterranean Sea and modelled up to a depth of 50 km. In this section, we first derive the accuracy of available gravity and magnetic models and then present the realistic synthetic model used to set up the inversion.
Starting from gravity data, we used the XGM2019e model (Zingerle et al., 2020) and to assess its accuracy we exploited different datasets: a freely available set of shipborne gravity data retrieved from the International Gravimetric Bureau (BGI), a compilation of gravity data on the Tyrrhenian and Adriatic regions from Zahorec et al. (2021) and a set of proprietary data (derived from shipborne surveys) available at Geomatics Research & Development srl (GReD). The accuracy of the gravity data from BGI and Zahorec et al. (2021) is provided by the institute that delivered the data, while for shipborne data available in GReD, having also raw observations, we inferred it from the analysis of crossovers following the procedure described in Sampietro et al. (2017). The XGM2019e model in the Mediterranean area has a standard deviation of about 3 mGal in terms of gravity anomalies, which gives about 7 E in terms of the second radial derivative (Murböck, 2015). Starting from the comparison with the available shipborne datasets, we also computed (see Figure 1) the error covariance of the XGM2019e model both in terms of gravity anomalies and second radial derivative (see, e.g., Hofmann-Wellenhof & Moritz, 2006, for the definition of the two functionals).
For both gravity functionals, a spatial correlation of the order of 15 km is found. Considering that the available proprietary shipborne data have an accuracy of about 0.5 mGal, we can assert that the XGM2019e accuracy estimate on the Mediterranean Sea area is reliable.
For the magnetic field, we started our analysis looking to available global scalar anomaly intensity grids such as the EMAG2v3 model (Meyer et al., 2016) and the WDMAM-2 (Lesur et al., 2016). Similarly to the comparison performed for the gravitational field, we compared EMAG2v3 and WDMAM-2 models with available proprietary magnetic data. The results are shown in terms of covariances in Figure 2.
In particular, the signal covariances of the two models on the Mediterranean Sea area are reported together with the covariances of the difference between EMAG2v3 and WDMAM-2 and between each global model and the available local data. The two global grids show differences of about 30 nT in terms of standard deviation (STD) over a signal of the order of 54 nT (STD). The comparison between the two global grids with the available local data shows an STD of the difference of about 23 nT. In the following, we will use the EMAG2v3 model assuming an error STD of 23 nT.

Realistic synthetic model
We then built the first synthetic model on the Mediterranean Sea area based on the following discontinuity surfaces: 1. bathymetry, 2. base of Plio-Quaternary sediments, 3. base of Messinian sediments, 4. basement and 5. Moho.
The bathymetry is taken from the GEBCO model (Mayer et al., 2018); the base of Plio-Quaternary and Messinian sediments has been derived by digitizing available models and seismic profiles (Capponi et al., 2020;de Voogd et al., 1992;Feld et al., 2017;Haq et al., 2020;Longacre et al., 2007;Makris & Yegorova, 2006). As for the basement, we derived it from the total sediment thickness of the World's Oceans and Marginal Seas (Version 3), developed by the National Oceanic and Atmospheric Administration (Straume et al., 2019), which, however, still largely depends on the Mediterranean Sea area from the CRUST2.0 global dataset (Molinari & Morelli, 2011 In the absence of better information, we assume the Curie isotherm depth to be equal to the Moho discontinuity. Characteristics for the layers are reported in Table 1, and the respective covariance functions are reported in Figure 3. The reference densities and magnetic susceptibilities are reported in Table 2. The crustal and upper mantle density models are taken from Fullea et al. (2021), while the susceptibility is derived from Hunt et al. (1995). In particular, crust susceptibility is obtained by scaling the densities from Fullea et al. (2021) to the value reported in Hunt et al. (1995). Starting from these reference sources, we built a threedimensional (3D) model with a spatial resolution of 1-arcmin (about 1.6 km) in and , 100 m in the direction up to a depth of 50 km for a total of about 1.14 × 10 9 cells. The expected accuracy of the shallower layers (base of Plio-Quaternary and base of Messinian sediments) is of few hundred metres while that of the basement and Moho is of few kilometres (as indicated in the literature used to build F I G U R E 3 Covariance functions of the reference geological horizons the model). Comparing these numbers with the STD of the amplitude of the different layers reported in Table 2, we can expect relative errors of about 20-30% for the shallowest layers and up to 50% for the deepest ones. A more difficult task is to quantify the accuracy of densities and susceptibilities. Actually, considering the range of variability at a global level of these parameters for different sedimentary and crustal rocks (see, e.g., Tenzer & Gladkikh, 2014;Christensen & Mooney, 1995;Hunt et al., 1995), we can expect errors in the model smaller than 5-10%, a part for the crustal susceptibility which strictly depends on the nature of the crust itself, and can also vary of the orders of magnitude. In any case, we want to emphasize that, since we will evaluate only different approximations of the above model, it is enough that the model is realistic, as it is, being derived from previously published studies.

Selection of observations ( , , )
We will start our discussion analysing the sensitivity of the different observables (namely gravity anomalies, second radial derivatives of gravitational potential and total magnetic intensity) to the main geological layers in order to choose the optimal one. In this study, we will consider for the gravity, the classical gravity anomalies ( ) and the second radial derivative of the potential ( rr ), both computed at sea level; for the magnetic field, the total magnetic intensity dF at 4 km above the surface. Figure 4 shows the computed gravitational effect of the synthetic model for the two gravity functionals (left panels) and the observed gravity field for the same functionals (right panels), while Figure 5 shows the effect of the synthetic model in terms of total magnetic intensity anomaly (left) and the observed total magnetic intensity anomaly (right). The observed gravity field has been synthesized from the XGM2019e model and reduced for the effect of the mantle (from 50 to 300 km) exploiting the WINTERC-G model (Fullea et al., 2021). As for the magnetic field, we used the EMAG2v3 model. A simple visual comparison with the actual observed fields shows that the computed forwards have the same order of magnitude and similar wavelengths of the gravity and magnetic observations. We then compute the gravitational effect of each geological horizon and compare the signal-to-noise ratio (Table 3) for rr , and dF. The signal-to-noise ratio, that is, the ratio between the signal STD and the observation error STD, can be used to understand the capability of a given functional to detect variations of a geological horizon. A signal-to-noise ratio smaller than 1 means that the observation error is larger than the signal and, therefore, that the specific functional cannot detect variations in the studied layer. Moreover, the larger is the signal-to-noise ratio the better the specific functional can detect variations on the geological horizon.
We see from Table 3 that gravity anomalies are better than the second radial derivative (having a higher signal-to-noise ratio) for every considered layer. Of course, this is true on the whole Mediterranean Sea area while for more local regions the situation could be different and the signal-to-noise ratio for rr could be larger than the one of .

F I G U R E 4 Computed gravitational effect of the synthetic model in terms of
(upper left panel) and rr (lower left panel) and observed gravity field in terms of (upper right panel) and rr (lower right panel)

F I G U R E 5 Computed total magnetic intensity anomaly due to the synthetic model (left panel), and observed total magnetic intensity (right panel)
Being interested in the whole Mediterranean Sea area, our analysis will be performed assuming as gravity observations. From Table 3, we also observe that, with the available global magnetic field data, we will be able to infer information only on the crystalline crust structure: the effect of this layer shows in fact a signal STD of about 37 nT. As for the other layers, only the Plio-Quaternary sediments, with the susceptibility considered in the synthetic model, locally reach a magnetic effect detectable by the accuracy of available global grids showing a signal-to-noise ratio larger than 1 close to the Nile Delta where the signal has its maximum.

Planar versus spherical approximation
We now rely on the realistic synthetic model, its gravity and magnetic predictions, and the true observation errors to investigate the effect of the planar approximation over this large region. We compute the gravitational effect of the synthetic model in planar and spherical approximations and compare the differences between the two with the corresponding observation error. The computation of the gravitational effect in the spherical approximation of a given volume is performed by means of tesseroids (Baykiev et al., 2016;Uieda et al., 2016). For the planar approximation, we use the closed solution of a right rectangular prism from Nagy (1966) and Bhattacharyya (1964). Figure 6 shows the difference between the forward responses of reference volume in planar and spherical approximations for and rr . Since, as stated in the second section , the average field does not enter in the inversion strategy, we omit it and plot the differences after the removal of their average values.
The signal is dominated by the border effect, clearly visible with typical bell-shaped anomalies. These effects are easily F I G U R E 6 Difference between the gravitational effect of the reference volume in planar and spherical approximations in terms of (left) and rr (right)

F I G U R E 7
Difference between the gravitational effect of the reference volume in planar and spherical approximations in terms of (left) and rr (right) once border effects are removed removed both in planar and spherical approximations by assuming the borders to have a constant densitȳ: • in the planar approximation, this is done by computing the effect of a Bouguer slab with densitȳand thickness equal to 50 km and removing the effect of a prism with the same side of our reference volume and densitȳ. • in the spherical approximation, the Bouguer slab is substituted by a spherical shell and the prism by a tesseroid.
Once removed the border effect (here we used as̄the average density of the synthetic model, that is 2938 kg/m 3 ), we obtain the results shown in Figure 7. We notice that the planar approximation in terms of entails errors of about 2 mGal (standard deviation) even considering the whole Mediterranean Sea area, which are smaller than the observation error. The error is largely correlated with the observed signal and is dominated by wavelengths larger than 1000 km. As a consequence in the case of the planar approximation when using only a small distortion of the deepest layers (and density distribution) is expected (e.g. by applying Parker-Oldenburg inversion (Oldenburg, 1974), we see that a Moho undulation with an STD of only 70 m can explain this signal). Moving to the second radial derivative instead, even after correcting the edge effects, the modelling error remains larger than the one of the observations with high-frequencies signals. The same procedure is also applied F I G U R E 8 Difference between the total magnetic intensity of the reference volume in planar and spherical approximations to estimate the planar approximation error for the magnetic field. The difference between planar and spherical approximations, supposing for the sake of simplicity an induced field directed in the radial and up directions, has a standard deviation of about 13 nT as shown in Figure 8. Similar results are also obtained for different directions of the inducing field as far as they are correctly projected between the spherical and planar reference systems.

Border size for inversion
An important point to be considered with potential field inversion is related to the proper dimension of the border around the volume of interest to be considered when performing the inversion. As shown in Figure 6, if not considered, edge effects can be larger (even of one order of magnitude) than the observed gravity anomaly. The size of the border is estimated empirically considering a target accuracy given by the observation error, that is, 3 mGal − 23 nT for our test case. We start increasing the size of the border, in such a way that the STD of the gravitational and magnetic fields above the volume (where we have the observations) is reduced. We select the minimum border giving an STD smaller than the target accuracy. An efficient way to compute the effect of the border is by means of a mixed solution (Sampietro & Capponi, 2019) where part of the borders is taken from global crustal models (such as CRUST1.0 Laske et al. (2013) or LITHO1.0 Pasyanos et al. (2014) or WINTERC- G Fullea et al. (2021)) and the remaining part by extending the reference density distribution values at the lateral boundaries of the volume by means of the nearest neighbour interpolation. This interpolation presents the big advantage of giving a density model of the borders which is continuous and allows to reduce the computational time required to evaluate its effect. In our simple test, we will use constant densities again by computing the effect with tesseroids for the spherical approximation and prism for the planar one. We show in Figure 9 the results obtained for the planar approximation in terms of , rr and dF. Starting from the gravity field, we observe that a border of about 15 • should be considered in order to reduce edge effects STD to values lower than the observation error. As expected, this value is notably reduced when moving to rr to less than 5 • . A behaviour similar to the one of rr is found for the magnetic field. However, this time, due to the high observation error STD, a border of about 1 • is sufficient to reduce edge F I G U R E 1 0 Error (mantle effect) as a function of the maximum depth in terms of effect STD to less than 23 nT. Note that, in case of good quality magnetic observation, with 1 nT accuracy a border of about 5 • is required. This is confirmed by a study on the required border at 4 km altitude by Baykiev et al. (2016) where the authors used magnetic forward in terms of tesseroids to compute the radial component of the crustal magnetic field on two areas for different border size. The situation is also very similar when moving to the spherical approximation, with being the functional requiring the largest border size. What we found is that, in general, for the spherical approximation, an increase of the size of the border up to more than 30 • for and 10 • for rr is required.

Maximum depth of the three-dimensional model
Similarly to the border size in the latitude, longitude (northsouth, west-east) directions, in order to deal with gravity, we also have to analyse the maximum required depth of the considered three-dimensional (3D) volume. This is not necessary when working with the magnetic field only because below the Curie isotherm no magnetization (except for the main field contribution) is expected. Supposing that the mantle density distribution is known with sufficient accuracy (see Sampietro & Capponi, 2019) and that we are interested in inverting the volume up to a depth of 50 km (slightly more than the expected Moho), we simply use a global mantle density model to strip this effect from the observations. In detail, the WINTERC-G model has been used, interpolating the density distribution of the mantle from a depth of 50 km down to 400 km with a discretization in the direction of 1 km. In Figure 10, the possible omission error (STD) from neglecting the mantle as a function of the maximum considered depth is shown. The plot shows, for instance, that removing from the data the effect of the density distribution between 50 and 100 km only entails an omission error of about 30 mGal. It turns out that, for the considered region, the effect of the mantle at least up to about 300 km should be removed from the observations. We also found that the effect of the density variation between

Model resolution
We also use our synthetic model to optimize the 3D model spatial resolution. Here we propose a new simple (empirical) approach where we compare the fields generated by the full 3D model with that obtained by downsampling the synthetic model of a factor . We increase as long as the downsampling error is smaller than the observed one. We report in Figure 11 maps of the downsampling error in terms of for = 2, 5, 20, 40 corresponding to a spatial resolution of 3.2, 8, 32, 64 km (recall that the initial model has been discretized in cells of size 1.6 km×1.6 km in the planar direction). Figure 12 shows the STD of this downsampling error in the and directions. In the figure, the observation error for and rr is also reported. We see how, considering the whole Mediterranean Sea area and relying on the XGM2019e and the EMAG2v3 models = 20, corresponding to a spatial resolution of 32 km seems to be appropriate.
In fact, for this spatial resolution, the total error is still dominated by the observation error and therefore basically no degradation of the initial information content occurred. Note also that a downsampling factor = 20 will reduce the total number of cells of our 3D model from about 1.14 × 10 9 to about 2.84 × 10 6 cells. When we move to higher values of , the total error starts to be dominated by the model one, thus F I G U R E 1 2 Error from downsampling in the and directions for the considered functionals. Dashed lines represent observation errors for (blue) and rr . Units are mGal, E, and nT for , rr and dF, respectively. not allowing to completely exploit the observation content. This is true whatever inversion algorithm is applied since in the end the total error will be interpreted by the inversion in terms of density distribution.
The same exercise is also performed to set the size of the volumetric elements in the vertical direction. Again, starting from the initial resolution of 100 m, we increased the size of volumetric elements and evaluated the difference between the initial and the simplified models. Results are reported in Figure 13, where it can be seen that a resolution of about 600 m is enough to discretize the volume in the vertical direction.

SYNTHETIC CASE STUDY
In the previous sections, we have shown how a wrong set of the considered parameters will lead to a distortion in the observations. This distortion in the observations is reflected, as obvious, as a distortion of the final results. We also provide empirical schemes to deal with this problem. To further substantiate the logic behind the proposed methodology, we performed a (non-exhaustive) set of inversions on a synthetic case study. In doing this, we have to consider that, in a potential field inversion (even on a synthetic model), possible errors in the results due to the mismodelling of the parameters considered in our work are mixed up with instability and nonuniqueness problems typical of the specific inversion scheme used. To remove at least the non-uniqueness, we keep the focus on the Mediterranean Sea area but simplify the volume to a two-layer problem, that is to the estimate of the Moho discontinuity supposing perfectly known all the other quantities. This simplification guarantees the uniqueness of the inverse problem solution (Sampietro & Sansò, 2012) and allows an easier evaluation of the results which can be assessed by comparing the true reference Moho with the estimated one. In this simplified model, we consider a standard density of 2670 kg/m 3 for the crust and of 3330 kg/m 3 for the mantle separated by the Moho model presented in the third section. The true model is created with a cell size of 1.6 km in the planar direction and 500 m in the vertical one and adding a border of 4000 km. A 3-mGal correlated noise, with a correlation length of 16 km has also been simulated. The gravitational effect due to the true model and the simulated observation error is reported in Figure 14.
As the inversion scheme, we used the Bayesian algorithm described in Sansó and Sampietro (2022). Since this algorithm is based on a Bayesian assumption, it requires an a priori model of the studied volume. We used as prior a simple model obtained by supposing a flat Moho at the average depth of the true Moho (i.e. about 31 km) and allowing to change this a priori surface between the maximum and the minimum depth of the true Moho. The true and the a priori models are shown in Figure 15.
We start by evaluating the effect of the border size. This is done by performing the inversion with borders of 4000, 1500, 1000 or 500 km. The difference between the true Moho and the one obtained by the inversion considering the 4000-km border size is shown in Figure 16 and has an STD of 396 m.
Since the true model has been built exactly with a border of 4000 km, the differences in Figure 16 are basically the combined effect of observation errors and of errors due to the specific inversion used. Figure 17 shows the differences between the true and estimated Moho for the other border sizes.
For borders smaller than about 1500 km, the solution is degraded with a worsening of the results of about 110% (for 500 km border) and 27% (for 1000 km border). Considering the 500-km border solution, the difference between the true and the estimated Moho is dominated by long wavelength effects. The larger is the considered border the smallest are these effects. When the border reaches the 1500-km, the residuals show the same order of magnitude of the 4000-km border solution. So has claimed in the paper larger borders would not significantly improve the inversion result.
In the second test, we consider the model spatial resolution in the planar direction. We perform different inversions with spatial resolutions of 8, 33 and 66 km. Looking to the differences with respect to the true model (Figure 18), we see that the 8 and 33-km spatial resolutions show almost the same error with STD of 394 and 402 m.
When moving to the 66-km spatial resolution, a degradation of the retrieved Moho of about 18% is obtained with an STD of the differences of more than 450 m. Main differences are observed in the Black Sea, close to Sicily or between Sardinia and Balearic Islands. We notice that in terms of cell numbers, we move from about 8 × 10 6 cells of the 8-km spatial resolution model to the 5 × 10 5 of the 33-km with a great saving of memory and computational time. Also in this case, results are well aligned with the spatial resolution suggested by the proposed methodologies which predicted negligible effects with about 30 km spatial resolution.
What we can conclude from these numerical tests, which as said are not exhaustive of the whole aspects considered in our work, is that the results obtained with the inversions on the simplified model are well aligned to those predicted by the proposed methodology, thus providing a further empirical validation to the logic behind our work.

CONCLUSIONS
The scope of the current work was to define an empirical procedure to calibrate all the ancillary parameters required before performing potential field inversion. The proposed methodology mainly relies on the comparison between the observation error (considered as a fixed threshold) and model errors, In the work, the Mediterranean Sea area has been used as a test case, with the dual objective to practically show how to calibrate the ancillary parameters and to present a benchmark for a future study dealing with regional areas. Looking to the results in the Mediterranean Sea region, we first assessed the accuracy of global gravity and magnetic field models: the XGM2019e global gravity field model in terms of synthesized at sea level shows an accuracy of 3 mGal, while the global magnetic grid from EMAG2v3 at 4 km altitude has an accuracy of 20 nT. After that, we built a realistic synthetic three-dimensional model in terms of geological units, densities and susceptibilities on the considered volume. Relying on the above observation errors and synthetic model, we performed a set of simulations to prove that even for such a large region the planar approximation is sufficient entailing errors smaller than the observations ones. Looking to the F I G U R E 1 8 Differences between the true and reference Moho for different spatial resolutions signal-to-noise ratio for the different geological horizons of the synthetic model, we also observed that, on the whole Mediterranean Sea area, for any horizon, gravity anomalies are the preferred functional to be used with respect to second radial derivatives (when both functionals are synthesized at ground level) because they have a higher ratio, meaning that they have more potential to detect variations within a certain geological horizon. Considering a maximum depth of the inversion volume of 50 km, borders of 15 • , roughly corresponding to about 1600 km are required to ensure a border effect below the observation error threshold. A model resolution in the and directions of about 30 km is sufficient, while for the vertical direction the suggested discretization is of the order of 600 m. The gravitational effect of the upper mantle should be modelled up to a depth of at least 300 km (e.g. from the WINTERC-G model) and removed from the observed field. From the computational point of view, the most demanding operation is the comparison between the planar and the spherical approximations. The latter requiring to compute the effect of a large and complex model by means of tesseroids. With the solution used in the paper, computed on a standard desktop PC (Intel(R) Core(TM) i7-4790K CPU 4.00GHz with 32 GB RAM), based on Gauss quadrature integration a computational time of about 8 h is required. However, the application of well-known acceleration by means of one-dimensional fast fourier transform (FFT) and the use of parallel computing are expected to reduce the computational time by more than a factor 30 (see, e.g., Zeng et al. 2022). Once the applicability of the planar approximation has been proven, the forward of our high-resolution model on the Mediterranean Sea area, performed by exploiting the right rectangular prism and FFT, requires less than 5 s, thus allowing to set all the parameters investigated in less than 10 min.
Some inversion tests, performed on a synthetic case study, give a further validation to the logic behind the proposed method, showing results well aligned with to those predicted by our empirical procedure to calibrate the inversion ancillary parameters.
The main limitation of the proposed procedure is that it needs the proper knowledge of the volume investigated and of the error of the available observations. The former is usually available at least roughly from the literature or in the worst case from global crustal models (such as the WINTERC-G model). The latter, if not known, can be obtained by comparing different datasets in the same area (as in the current work), or can be estimated starting from a literature review and knowing how the observations have been acquired and processed.

A C K N O W L E D G E M E N T S
This study was supported by European Space Agency, ESA Contract No. 400013688421I-DT-lr.

D A T A AVA I L A B I L I T Y S T A T E M E N T
The main results of this study, carried out in the framework of the European Space Agency XORN project, are freely available, upon request to the authors O R C I D Daniele Sampietro https://orcid.org/0000-0001-9747-8497 Erwan Thébault https://orcid.org/0000-0002-7038-2855 Lydie Gailler https://orcid.org/0000-0002-8132-2428