<figure><img src="image_1.png" /></figure>

HAL
open science

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

Daniele Sampietro, Martina Capponi, Erwan Thébault, Lydie Gailler

► To cite this version:

Daniele Sampietro, Martina Capponi, Erwan Thébault, Lydie Gailler. An Empirical Method for the Optimal Setting of the Potential Fields Inverse Problem. Geophysical Prospecting, 2023, 71 (2), pp.350-365. 10.1111/1365-2478.13305. hal-03872970

HAL Id: hal-03872970

https://hal.science/hal-03872970

Submitted on 26 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
An Empirical Method for the Optimal Setting of the Potential
Fields Inverse Problem

<table><tr><td>Journal:</td><td>Geophysical Prospecting</td></tr><tr><td>Manuscript ID</td><td>GP-2022-3988.R4</td></tr><tr><td>Wiley - Manuscript type:</td><td>Original Article</td></tr><tr><td>Date Submitted by the Author:</td><td>21-Nov-2022</td></tr><tr><td>Complete List of Authors:</td><td>Sampietro, Daniele; Geomatics Research &amp; Development srl, Geomatics Research &amp; Development srl<br>Capponi, Martina; Geomatics Research &amp; Development srl, Geomatics Research &amp; Development srl<br>Thébault, Erwan; Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans<br>Gailler, Lydie; Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans</td></tr><tr><td>Keyword:</td><td>Potential field inversion, 3D problem, Gravity modeling, Magnetic data analysis</td></tr></table>
# An Empirical Method for the Optimal Setting of the Potential Fields Inverse Problem

Daniele Sampietro¹ | Martina Capponi¹ | Erwan
Thébault² | Lydie Gailler²

¹Geomatics Research & Development srl, Lomazzo (CO), Italy

²Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, Clermont-Ferrand, France

**Correspondence**

Sampietro D., Geomatics Research &
Development srl, Lomazzo (CO), 22074,
Italy

Email: daniele.sampietro@g-red.eu

**Funding information**

European Space Agency, ESA Contract No.
400013688421I-DT-Ir

The use of potential field methods for geophysical explo-
ration purposes is nowadays quite common: these tech-
niques consent to retrieve geological knowledge over ex-
tended regions and can give complementary information
where other invasive or expensive techniques, such as seis-
mic 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 sur-
veys, 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, i.e. about 3 – 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 regional scale. In the
present work some preliminary technical aspects related to
the use of these models to perform 3D inversion are dis-
cussed, thus defining an empirical but rigorous procedure
to setup gravity and magnetic inversion. In particular, we
address the questions whether the classical planar approx-
imation is acceptable for regional inversions or if a spheri-
cal 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 3D volume area to be modeled depending on the specific target investigated. The application of the proposed methods to the Mediterranean case study is also presented.

KEYWORDS

Gravity and Magnetic Inversion, 3D inversion, parameters setting, Mediterranean Sea

# 1 | 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 near surface regional scale. Inverse gravimetric/magnetic problems aim at reconstructing subsurface properties given the observed field above the Earth 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 suggesting the papers of Nabighian et al. (2005a,b) 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. Li and Oldenburg (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 e.g. 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 a 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;

* 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. 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. 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 (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, 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 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 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.

## 2 | 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 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. Li and Oldenburg (1996)):

$$ \phi = [y_o - F(x)]^T C_{vv}^{-1} [y_o - F(x)] \quad (1) $$

where $y_o$ is the vector of the observations, $x$ is the vector of unknowns, $C_{vv}$ 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):

$$ F_{\delta g} = -G \iiint \rho(x, y, z) \frac{z - z_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2}} dx dy dz \quad (2) $$

where $x_i$, $y_i$, $z_i$ are the coordinates of the i-element of the observations vector, $G$ is the Newton's gravitational constant and $x$, $y$ and $z$ 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):

$$ F_B = -C_m B \chi(x, y, z) \cdot \nabla_P \iiint \frac{B \chi(x, y, z) \cdot \nabla_Q}{\sqrt{(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2}} dx dy dz \quad (3) $$
where C<sub>m</sub> is a constant depending on the used unit system, B 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 V should be limited both in the lateral extent and in depth. Moreover ρ(x, y, z) and χ(x, y, z) 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 in dividing the whole volume V in a set of rectangular prisms with size Δx × Δy × Δz.

Considering as observations in Eq. 1 the gravity and magnetic fields, and assuming measurement errors for the two quantities independent with covariances C<sub>δgδg</sub> and C<sub>BB</sub> respectively, Eq. 1 becomes:

$$ \phi = [y_{\delta g} - \bar{F}_{\delta g} \rho]^T C_{\delta g \delta g}^{-1} [y_{\delta g} - \bar{F}_{\delta g} \rho] + [y_B - \bar{F}_B \chi]^T C_{BB}^{-1} [y_B - \bar{F}_B \chi] \quad (4) $$

where $\bar{F}_{\delta g}$ is a $m_{\delta g} \times n$ matrix, with $m_{\delta g}$ number of gravity observations and $n$ number of discretized densities, containing at the element $i, j$ the gravitational effect of the prism $j$ to the observation point $i$ (obtained by integrating $F_{\delta g}$). Similarly, $\bar{F}_B$ is a $m_B \times n$ matrix, obtained by integrating $F_B$. Eq. 4 represents the base of every joint inverse problem. At this point, we are left with a set of given observations $y_{\delta g}$, $y_B$ to be inverted in order to estimate a finite set of densities ρ and magnetic susceptibilities χ, within the discetized study region V. In such a problem we can consider as known the observations as well as their accuracy and spatial resolution. On the contrary, the volume V (defined by the coordinates of its vertices) and the size of the prisms in which it is discretized (namely Δx, Δy and Δz) should be defined before the inversion. In addition, to avoid edge effects, an external border region of side Δb should be also considered.

An important remark has to be made on the forward operators defined in Eq. 4. When dealing with such regional inversion problems, the gravitational effect of the Earth 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 and Capponi, 2019). Choosing the latter possibility, the first term of the data-misfit equation can be written as:

$$ \phi_{\delta g} = [y_{\delta g} - E[y_{\delta g}] - \left(I - \frac{J}{m_{\delta g}}\right) \bar{F}_{\delta g} \rho]^T C_{\delta g \delta g}^{-1} [y_{\delta g} - E[y_{\delta g}] - \left(I - \frac{J}{m}\right) \bar{F}_{\delta g} \rho] \quad (5) $$

where E[y<sub>δg</sub>] is the expectation operator applied to y<sub>δg</sub>, $m_{\delta g}$ is the number of observations, I is the identity matrix of size $m_{\delta g}$ and J is an $m_{\delta g} \times m_{\delta g}$ all-ones matrix. Basically $(I - \frac{J}{m_{\delta g}}) \bar{F}_{\delta g}$ represents a zero-mean forward operator (Sansó and Sampietro, 2022).

A similar argument can be used also for the magnetic field (at least for regions smaller than 2500 km × 2500 km). In fact, when trying to recover information on the Earth crust from the magnetic field, the lithospheric contribute 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., 2016, 2021). Total magnetic internal field models show a break in the slope of degree variances between degrees 13 and 16 (Langel and Estes, 1982) roughly corresponding to wavelengths between 3000 km 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 happen 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 implies 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 parametrization of the inverse problem, we will be lead 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 affect 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 fields 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 $(I - \frac{J}{m\delta_g})\bar{F}_{\delta_g}$ and $(I - \frac{J}{m_B})\bar{F}_B$ operators. To evaluate these model errors, the effect of the studied volume (in terms of the grav/mag potential field) will be firstly 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.

## 3 | 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 firstly derive the accuracy of available gravity and magnetic models and then present the realistic synthetic model used to setup 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 who delivered the data, while for shipborne data available in GReD, having also raw observations, we inferred it from the analysis of crossovers following to 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 of second radial derivative (see e.g. Hofmann-Wellenhof and 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
<figure><img src="image_3.png" /><figcaption>FIGURE 1 XGM2019e error covariance in the Mediterranean Sea area</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2 EMAG2v3, WDMAM-2 signal and error covariances.</figcaption></figure>

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

## 3.1 | Realistic Synthetic model

We then built a first synthetic model on the Mediterranean Sea area based on the following discontinuity surfaces:

1. bathymetry;

2. base of Plio-quaternary sediments;
<figure><img src="image_5.png" /><figcaption>FIGURE 3 Covariance functions of the reference geological horizons</figcaption></figure>

<table><caption>TABLE 1 Statistics of the reference geological horizons.</caption><thead><tr><th rowspan="2">LAYER</th><th>AVERAGE DEPTH</th><th>DEPTH STD</th><th>CORRELATION LENGTH</th></tr><tr><th>[m]</th><th>[m]</th><th>[m]</th></tr></thead><tbody><tr><td>Bathymetry</td><td>645</td><td>1090</td><td>15400</td></tr><tr><td>Base Plio-Quaternary</td><td>845</td><td>1274</td><td>18700</td></tr><tr><td>Base Mess. Salt</td><td>947</td><td>1647</td><td>17600</td></tr><tr><td>Basement</td><td>3943</td><td>3790</td><td>28600</td></tr><tr><td>Moho</td><td>31218</td><td>6810</td><td>33000</td></tr></tbody></table>

3. base of Messinian sediments;

4. basement;

5. Moho;

The bathymetry is taken from the GEBCO model (Mayer et al., 2018), the base of Plio-Quaternary and Messinian sediments have been derived by digitizing available models and seismic profiles (Capponi et al., 2020; Haq et al., 2020; de Voogd et al., 1992; Longacre et al., 2007; Feld et al., 2017; Makris and 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 NOAA (Straume et al., 2019), which however still largely depends in the Mediterranean Sea area from the CRUST2.0 global dataset (Molinari and Morelli, 2011). Finally we took the Moho depth from the Moho Map of the European Plate (Grad et al., 2009). In absence of better information we assume the Curie isotherm depth to be equal to the Moho discontinuity. Characteristics for the layers are reported in Tab. 1 and the respective covariance functions are reported in Fig. 3. The reference densities and magnetic susceptibilities are reported in Tab. 2. The crustal and upper mantle density models are taken from Fullea et al. (2021), while the susceptibility is derived from Hunt et al.
<table><caption>**TABLE 2** reference density and susceptibility model. * Crust and Upper mantle models lateral density and susceptibility distribution have been taken from Fullea et al. (2021) and Hunt et al. (1995).</caption><thead><tr><th>LAYER</th><th>AVERAGE DENSITY<br>[kg/m³]</th><th>DENSITY GRADIENT<br>[kg/m³/km]</th><th>SUSCEPTIBILITY<br>[10<sup>-6</sup>SI]</th></tr></thead><tbody><tr><td>Water</td><td>1030</td><td>0</td><td>-13</td></tr><tr><td>Plio-Quaternary</td><td>2220</td><td>0</td><td>10000</td></tr><tr><td>Messinian</td><td>2170</td><td>0</td><td>-30</td></tr><tr><td>Pre-Messinian</td><td>2400</td><td>6</td><td>1500</td></tr><tr><td>Crust</td><td>2890*</td><td>6</td><td>36796*</td></tr><tr><td>Upper Mantle</td><td>3321*</td><td>0</td><td>0</td></tr></tbody></table>

(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 3D model with spatial resolution of 1-arcminute (about 1.6 km) in x and y, 100 m in the z direction up to a depth of 50 km for a total of about 1.14 · 10⁹ cells. The expected accuracy of the shallower layers (base of Plio-quaternary and base of Messinian sediments) is of few hundreds metres while that of basement and moho is of few kilometres (as indicated in the literature used to build the model). Comparing these numbers with the STD of the amplitude of the different layers reported in Tab. 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 global level of these parameters for different sedimentary and crustal rocks (see e.g. Tenzer and Gladkikh (2014); Christensen and 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 vary also of order 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 previous published studies.

## 4 | RESULTS

### 4.1 | Selection of observations ($δg$, $T_{rr}$, dB)

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 ($δg$) and the second radial derivative of the potential ($T_{rr}$), both computed at sea level; for the magnetic field, the total magnetic intensity dB at 4 km above the surface. Fig. 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 Fig. 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 km 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
<figure><img src="image_6.png" /></figure>

**FIGURE 4** Computed gravitational effect of the synthetic model in terms of $\delta g$ (upper left panel) and $T_{rr}$ (lower left panel) and observed gravity field in terms of $\delta g$ (upper right panel) and $T_{rr}$ (lower right panel)

<figure><img src="image_7.png" /></figure>

**FIGURE 5** Computed total magnetic intensity anomaly due to the synthetic model (left panel), and observed total magnetic intensity (right panel).
<table><caption>TABLE 3 Signal-to-noise ratio for the different observables in the whole Mediterranean Sea area</caption><thead><tr><th>LAYER</th><th>δg</th><th>Trr</th><th>dF</th></tr></thead><tbody><tr><td>Water</td><td>29</td><td>4</td><td>2·10<sup>-4</sup></td></tr><tr><td>Plio-Quaternary</td><td>10</td><td>4</td><td>0.1</td></tr><tr><td>Messinian</td><td>14</td><td>2</td><td>3·10<sup>-4</sup></td></tr><tr><td>Pre-Messinian</td><td>90</td><td>8</td><td>0.04</td></tr><tr><td>Crust</td><td>252</td><td>11</td><td>2</td></tr><tr><td>Upper Mantle</td><td>203</td><td>5</td><td>0</td></tr></tbody></table>

that the computed forwards have 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 $T_{rr}$, $\delta g$ and $dF$. The signal-to-noise ratio, i.e. 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 second radial derivative (having an 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 $T_{rr}$ could be larger than the one of $\delta g$. Being interested in the whole Mediterranean Sea area our analysis will be performed assuming $\delta g$ as gravity observations. From Tab. 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 reaches 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.

## 4.2 | Planar vs. 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 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 spherical approximation of a given volume is performed by means of tesseract (Uieda et al., 2016; Baykiev 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 approximation for $\delta g$ and $T_{rr}$. Since, as stated in Sec. 2, 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 border effect, clearly visible with typical bell-shaped anomalies. These effects are
<figure><img src="image_8.png" /><figcaption>FIGURE 6 Difference between the gravitational effect of the reference volume in planar and spherical approximations in terms of $\delta g$ (left) and $T_{rr}$ (right).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 7 Difference between the gravitational effect of the reference volume in planar and spherical approximations in terms of $\delta g$ (left) and $T_{rr}$ (right) once border effects are removed.</figcaption></figure>

easily removed both in planar and spherical approximation by assuming the borders to have a constant density $\bar{\rho}$:

* in planar approximation, this is done by computing the effect of a Bouguer slab with density $\bar{\rho}$ and thickness equal to 50 km and removing the effect of a prism with the same side of our reference volume and density $\bar{\rho}$.

* in spherical approximation the Bouguer slab is substituted by a spherical shell and the prism by a tesseract.

Once removed the border effect (here we used as $\bar{\rho}$ the average density of the synthetic model, i.e. 2938 kg/m³), we obtain the results shown in Fig. 7. We notice that the planar approximation in terms of $\delta g$ 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 case of planar approximation when using $\delta g$ 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 applied also to estimate the planar approximation error for the magnetic field. The difference between planar and spherical approximation, 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 Fig. 8.Similar results are obtained also for different directions of the inducing field as far as they are correctly projected between the spherical and planar reference systems.
<figure><img src="image_10.png" /><figcaption>**FIGURE 8** Difference between the Total magnetic intensity of the reference volume in planar and spherical approximations.</figcaption></figure>

## 4.3 | 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 Fig. 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, i.e. 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 and 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 a 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 tesseracts for the spherical approximation and prism for the planar one. We show in Fig. 9 the results obtained for the planar approximation in terms of δg, $T_{rr}$ and dF. Starting from the gravity field, we observe that a border of about 15° should be consider in order to reduce edge effects STD to values lower than the observation error. As expected this value is notably reduced when moving to $T_{rr}$ to less than 5°. A behaviour similar to the one of $T_{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 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 tesseracts to compute the radial component of the crustal magnetic field on two areas for different border size. The situation is very similar also when moving to the spherical approximation, with δg 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 δg and 10° for $T_{rr}$ is required.
<figure><img src="image_11.png" /><figcaption>**FIGURE 9** Error (border effect) as a function of the size of the border for $\delta g$ (upper panel), $T_{rr}$ (middle panel) and $dF$ (lower panel).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIGURE 10 Error (mantle effect) as a function of the maximum depth in terms of δg.</figcaption></figure>

## 4.4 | Maximum depth of the 3D model

Similarly to the border size in the latitude, longitude (North-South, West-East) directions, in order to deal with gravity we also have to analyse the maximum required depth of the considered 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 contribute) is expected. Supposing that the mantle density distribution is known with sufficient accuracy (see Sampietro and 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 details 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 *z* direction of 1 km. In Fig. 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 km 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 399 km and 400 km has an STD of about 0.01 mGal, thus making the effect of the deeper density distribution negligible.

## 4.5 | Model resolution

We use our synthetic model also 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 *n*. We increase *n* as long as the downsampling error is smaller than the observation one. We report in Fig. 11 maps of the downsampling error in terms of δg for *n* = 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). Fig. 12 shows the STD of this downsampling error in the *x* and *y* directions. In the figure the observation error for δg and $T_{rr}$ is also reported. We see how, considering the whole Mediterranean Sea area and relying on the XGM2019e and the EMAG2v3 models *n* = 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 *n* = 20 will reduce the total number of cells of our 3D model from about 1.14·10⁹ to about 2.84·10⁶ cells. When we
<figure><img src="image_13.png" /><figcaption>FIGURE 11 Downsampling error in terms of $\delta g$ for $n = 2, 5, 20, 40$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 12 Error from downsampling in the x and y directions for the considered functionals. Dashed lines represent Observation errors for $\delta g$ (blue) and $T_{rr}$. Units are mGal, E, nT for $\delta g$, $T_{rr}$ and dF respectively</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**FIGURE 13** Error due to the downsampling in the *z* direction for the considered functionals. Dashed lines represent observation errors for *δg* and *Trr*. Units are mGal, E, nT for *δg*, *Trr* and dF respectively</figcaption></figure>

move to higher values of *n* the total error start to be dominated by the model one thus 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 performed also 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 model. Results are reported in Fig. 13 where it can be seen that a resolution of about 600 m is enough to discretize the volume in the vertical direction.

5 | 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, in 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 non-uniqueness 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 simplifying the volume to a two-layer problem, i.e. 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 and 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³ for the crust and of 3330 kg/m³ for the mantle separated by the Moho model presented in the Sec. 3. 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 been also simulated. The gravitational effect due to the true model and the simulated observation error are reported in Fig. 14. As 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
<figure><img src="image_16.png" /></figure>

**FIGURE 14** gravitational effect due to the mass distribution of the true model (left panel) and simulated observation error (right panel).

<figure><img src="image_17.png" /></figure>

**FIGURE 15** True model (up), and a-priori information for the inversion (down). Green lines represent the range of variability of the a-priori Moho.

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 Fig. 15.

We start by evaluating the effect of the border size. This is done by performing the inversion with borders of 4000 km, 1500 km, 1000 km 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 Fig 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 Fig. 16 are basically the combined effect of observation errors and of errors due to the specific inversion used. Fig. 17 shows the differences between the true and estimated Moho for the other border sizes. For border 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 wavelengths 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
<figure><img src="image_18.png" /><figcaption>FIGURE 16 Difference between the true and estimated Moho with 4000 km border.</figcaption></figure>

improve the inversion result.

In a second test we consider the model spatial resolution in the planar direction. We perform different inversions with spatial resolutions of 8 km, 33 km and 66 km. Looking to the differences with respect to the true model (Fig. 18), we see that the 8 km and 33 km spatial resolutions show almost the same error with STD of 394 m and 402. 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 concentrated in the Black Sea, close to Sicily or between Sardinia and Balearic Islands. We notice that in terms of cells number we move from about 8 · 10⁶ cells of the 8 km spatial resolution model to the 5 · 10⁵ 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.

6 | 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, computed with a set of forwards performed on a realistic synthetic model of the considered volume.

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 future study dealing with regional areas. Looking to the results in the Mediterranean Sea region we firstly assessed the accuracy of global gravity and magnetic field models: the XGM2019e global gravity field model in terms of δg synthesised 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 3D model in terms of geological units, densities and susceptibilities on the considered volume. Relying on the above observations errors and synthetic model we performed a set of simulation to prove that even for such a large region the planar approximation is sufficient entailing errors smaller than the observations ones. Looking to the signal-to-noise ratio for the different geological horizons of the synthetic model we also observed that, on the
<figure><img src="image_19.png" /><figcaption>**FIGURE 17** Differences between the true and reference Moho for different border sizes.</figcaption></figure>
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**FIGURE 18** Differences between the true and reference Moho for different spatial resolutions.</figcaption></figure>
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 x and y 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 tesseract. 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 hours is required. However the application of well known acceleration by means of 1D FFT and the use of parallel computing are expected to reduce the computational time of 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 right rectangular prism and FFT requires less than 5 seconds, thus allowing to set all the parameters investigated in less than 10 minutes.

Some inversions tests, performed on a synthetic case study, gives 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 a proper knowledge of the volume investigated and of the error of the available observations. The former is usually available at least in a rough way from 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 literature review and knowing how the observations have been acquired and processed.

## Acknowledgements

Acknowledgements should include contributions from anyone who does not meet the criteria for authorship (for example, to recognize contributions from people who provided technical help, collation of data, writing assistance, acquisition of funding, or a department chairperson who provided general support), as well as any funding or other support information.

## Conflict of Interest

Daniele Sampietro, Martina Capponi, Erwan Thebauld, Lydie Gailler declare that they have no conflict of interest.

## Data Availability Statement

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
references

An, M. (2012) A simple method for determining the spatial resolution of a general inverse problem. *Geophysical Journal International*, **191**, 849-864.

Backus, G., George, B., Parker, R. L., Parker, R. and Constable, C. (1996) *Foundations of geomagnetism*. Cambridge University Press.

Baykiev, E., Ebbing, J., Brönner, M. and Fabian, K. (2016) Forward modeling magnetic fields of induced and remanent magnetization in the lithosphere using tesseract. *Computers & geosciences*, **96**, 124-135.

Bhattacharyya, B. (1964) Magnetic anomalies due to prism-shaped bodies with arbitrary polarization. *Geophysics*, **29**, 517-531.

Blakely, R. J. (1996) *Potential theory in gravity and magnetic applications*. Cambridge university press.

Bosch, M., Meza, R., Jiménez, R. and Hönig, A. (2006) Joint gravity and magnetic inversion in 3d using monte carlo methods. *Geophysics*, **71**, G153-G156.

Capponi, M., GlaviCh, E. and Sampietro, D. (2020) Map of plio-quaternary sediment depths in the mediterranean sea. *Bollettino di Geofisica Teorica e Applicata*, **61**, 421-432.

Christensen, N. I. and Mooney, W. D. (1995) Seismic velocity structure and composition of the continental crust: A global view. *Journal of Geophysical Research: Solid Earth*, **100**, 9761-9788.

Feld, C., Mechie, J., Hübscher, C., Hall, J., Nicolaides, S., Gurbuz, C., Bauer, K., Louden, K. and Weber, M. (2017) Crustal structure of the eratosthenes seamount, cyprus and s. turkey from an amphibian wide-angle seismic profile. *Tectonophysics*, **700**, 32-59.

Fullea, J., Lebedev, S., Martinec, Z. and Celli, N. (2021) Winterc-g: mapping the upper mantle thermochemical heterogeneity from coupled geophysical-petrological inversion of seismic waveforms, heat flow, surface elevation and gravity satellite data. *Geophysical Journal International*, **226**, 146-191.

Galley, C. G., Lelièvre, P. G. and Farquharson, C. G. (2020) Geophysical inversion for 3d contact surface geometry. *Geophysics*, **85**, K27-K45.

Grad, M., Tiira, T. and Group, E. W. (2009) The moho depth map of the european plate. *Geophysical Journal International*, **176**, 279-292.

Haq, B., Gorini, C., Baur, J., Moneron, J. and Rubino, J.-L. (2020) Deep mediterranean's messinian evaporite giant: How much salt? *Global and Planetary Change*, **184**, 103052.

Hofmann-Wellenhof, B. and Moritz, H. (2006) *Physical geodesy*. Springer Science & Business Media.

Hunt, C. P., Moskowitz, B. M., Banerjee, S. K. et al. (1995) Magnetic properties of rocks and minerals. *Rock physics and phase relations: A handbook of physical constants*, **3**, 189-204.

Kuhn, M., Featherstone, W. and Kirby, J. (2009) Complete spherical bouguer gravity anomalies over australia. *Australian journal of Earth sciences*, **56**, 213-223.

Langel, R. and Estes, R. (1982) A geomagnetic field spectrum. *Geophysical Research Letters*, **9**, 250-253.

Laske, G., Masters, G., Ma, Z. and Pasyanos, M. (2013) Update on crust1.0—a 1-degree global model of earth's crust. In *Geophys. Res. Abstr*, vol. 15, 2658.

Lesur, V., Hamoudi, M., Choi, Y., Dyment, J. and Thébault, E. (2016) Building the second version of the world digital magnetic anomaly map (wdmam). *Earth, Planets and Space*, **68**, 1-13.
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

Li, X. (2001) Vertical resolution: Gravity versus vertical gravity gradient. *The Leading Edge*, **20**, 901–904.

Li, Y. and Oldenburg, D. W. (1996) 3-d inversion of magnetic data. *Geophysics*, **61**, 394–408.

– (1998a) 3-d inversion of gravity data. *Geophysics*, **63**, 109–119.

– (1998b) Separation of regional and residual magnetic field data. *Geophysics*, **63**, 431–439.

Longacre, M., Bentham, P., Hanbal, I., Cotton, J. and Edwards, R. (2007) New crustal structure of the eastern mediterranean basin: detailed integration and modeling of gravity, magnetic, seismic refraction, and seismic reflection data. In *EGM 2007 International Workshop: Innovation in EM, Grav and Mag Methods: A New Perspective for Exploration*, vol. 15, 18.

Makris, J. and Yegorova, T. (2006) A 3-d density–velocity model between the cretan sea and libya. *Tectonophysics*, **417**, 201–220.

Mayer, L., Jakobsson, M., Allen, G., Dorschel, B., Falconer, R., Ferrini, V., Lamarche, G., Snaith, H. and Weatherall, P. (2018) The nippon foundation–gebco seabed 2030 project: The quest to see the world’s oceans completely mapped by 2030. *Geosciences*, **8**, 63.

Meyer, B., Saltus, R. W. and Chulliat, A. (2016) Emag2-v3: A new global compilation of lithospheric magnetic anomalies. In *AGU Fall Meeting Abstracts*, vol. 2016, GP43A–1210.

Molinari, I. and Morelli, A. (2011) Epcrust: a reference crustal model for the european plate. *Geophysical Journal International*, **185**, 352–364.

Moritz, H. (1990) The figure of the earth: theoretical geodesy and the earth’s interior. *Karlsruhe: Wichmann*.

Murböck, M. (2015) *Virtual constellations of next generation gravity missions*. Ph.D. thesis, Technische Universität München.

Nabighian, M. N., Grauch, V., Hansen, R., LaFehr, T., Li, Y., Peirce, J. W., Phillips, J. D. and Ruder, M. (2005a) The historical development of the magnetic method in exploration. *Geophysics*, **70**, 33ND–61ND.

– (2005b) The historical development of the magnetic method in exploration. *Geophysics*, **70**, 33ND–61ND.

Nagy, D. (1966) The gravitational attraction of a right rectangular prism. *Geophysics*, **31**, 362–371.

Oldenburg, D. W. (1974) The inversion and interpretation of gravity anomalies. *Geophysics*, **39**, 526–536.

Pasyanos, M. E., Masters, T. G., Laske, G. and Ma, Z. (2014) Litho1. 0: An updated crust and lithospheric model of the earth. *Journal of Geophysical Research: Solid Earth*, **119**, 2153–2173.

Sampietro, D. (2015) Geological units and moho depth determination in the western balkans exploiting goce data. *Geophysical Journal International*, **202**, 1054–1063.

Sampietro, D. and Capponi, M. (2019) Practical tips for 3d regional gravity inversion. *Geosciences*, **9**, 351.

Sampietro, D., Capponi, M., Mansi, A., Gatti, A., Marchetti, P. and Sansò, F. (2017) Space-wise approach for airborne gravity data modelling. *Journal of Geodesy*, **91**, 535–545.

Sampietro, D. and Sansò, F. (2012) Uniqueness theorems for inverse gravimetric problems. In *VII Hotine-Marussi Symposium on Mathematical Geodesy*, 111–115. Springer.

Sansó, F. and Sampietro, D. (2022) Analysis of the gravity field.

Straume, E. O., Gaina, C., Medvedev, S., Hochmuth, K., Gohl, K., Whittaker, J. M., Abdul Fattah, R., Doornenbal, J. C. and Hopper, J. R. (2019) Globsed: Updated total sediment thickness in the world’s oceans. *Geochemistry, Geophysics, Geosystems*, **20**, 1756–1772.
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

Tenzer, R. and Gladkikh, V. (2014) Assessment of density variations of marine sediments with ocean and sediment depths. *The Scientific World Journal*, **2014**.

Thébault, E., Hulot, G., Langlais, B. and Vigneron, P. (2021) A spherical harmonic model of earth's lithospheric magnetic field up to degree 1050. *Geophysical Research Letters*, **48**, e2021GL095147.

Thébault, E., Vigneron, P., Langlais, B. and Hulot, G. (2016) A swarm lithospheric magnetic field model to sh degree 80. *Earth, Planets and Space*, **68**, 1-13.

Uieda, L., Barbosa, V. C. and Braitenberg (2016) Tesseracts: Forward-modeling gravitational fields in spherical coordinates. *Geophysics*, **81**, F41-F48.

de Voogd, B., Truffert, C., Chamot-Rooke, N., Huchon, P., Lallemant, S. and Le Pichon, X. (1992) Two-ship deep seismic soundings in the basins of the eastern mediterranean sea (pasiphae cruise). *Geophysical Journal International*, **109**, 536-552.

Zahorec, P., Papčo, J., Pašteka, R., Bielik, M., Bonvalot, S., Braitenberg, C., Ebbing, J., Gabriel, G., Gosar, A., Grand, A. et al. (2021) The first pan-alpine surface-gravity database, a modern compilation that crosses frontiers. *Earth System Science Data*, **13**, 2165-2209.

Zeng, X., Wan, X., Lin, M. and Wang, W. (2022) Gravity field forward modelling using tesseracts accelerated by taylor series expansion and symmetry relations. *Geophysical Journal International*, **230**, 1565-1584.

Zingerle, P., Pail, R., Gruber, T. and Oikonomidou, X. (2020) The combined global gravity field model xgm2019e. *Journal of Geodesy*, **94**, 1-12.
# An Empirical Method for the Optimal Setting of the Potential Fields Inverse Problem

Daniele Sampietro¹ | Martina Capponi¹ | Erwan
Thébault² | Lydie Gailler²

¹Geomatics Research & Development srl, Lomazzo (CO), Italy

²Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, Clermont-Ferrand, France

**Correspondence**

Sampietro D., Geomatics Research &
Development srl, Lomazzo (CO), 22074,
Italy

Email: daniele.sampietro@g-red.eu

**Funding information**

European Space Agency, ESA Contract No.
400013688421I-DT-Ir

The use of potential field methods for geophysical explo-
ration purposes is nowadays quite common: these tech-
niques consent to retrieve geological knowledge over ex-
tended regions and can give complementary information
where other invasive or expensive techniques, such as seis-
mic 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 sur-
veys, 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, i.e. about 3 – 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 regional scale. In the
present work some preliminary technical aspects related to
the use of these models to perform 3D inversion are dis-
cussed, thus defining an empirical but rigorous procedure
to setup gravity and magnetic inversion. In particular, we
address the questions whether the classical planar approx-
imation is acceptable for regional inversions or if a spheri-
cal 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 3D volume area to be modeled depending on the specific target investigated. The application of the proposed methods to the Mediterranean case study is also presented.

KEYWORDS

Gravity and Magnetic Inversion, 3D inversion, parameters setting, Mediterranean Sea

# 1 | 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 near surface regional scale. Inverse gravimetric/magnetic problems aim at reconstructing subsurface properties given the observed field above the Earth 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 suggesting the papers of Nabighian et al. (2005a,b) 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. Li and Oldenburg (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 e.g. 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 a 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;

* 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. 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. Literature
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

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 (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, 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 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 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.

## 2 | 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 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. Li and Oldenburg (1996)):

$$ \phi = [y_o - F(x)]^T C_{vv}^{-1} [y_o - F(x)] \quad (1) $$

where $y_o$ is the vector of the observations, $x$ is the vector of unknowns, $C_{vv}$ 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):

$$ F_{\delta g} = -G \iiint \rho(x, y, z) \frac{z - z_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2}} dx dy dz \quad (2) $$

where $x_i, y_i, z_i$ are the coordinates of the $i$-element of the observations vector, $G$ is the Newton's gravitational constant and $x, y$ and $z$ 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):

$$ F_B = -C_m B \chi(x, y, z) \cdot \nabla_P \iiint \frac{B \chi(x, y, z) \cdot \nabla_Q}{\sqrt{(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2}} dx dy dz \quad (3) $$
where C<sub>m</sub> is a constant depending on the used unit system, B 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 V should be limited both in the lateral extent and in depth. Moreover ρ(x, y, z) and χ(x, y, z) 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 in dividing the whole volume V in a set of rectangular prisms with size Δx × Δy × Δz.

Considering as observations in Eq. 1 the gravity and magnetic fields, and assuming measurement errors for the two quantities independent with covariances C<sub>δgδg</sub> and C<sub>BB</sub> respectively, Eq. 1 becomes:

$$ \phi = [y_{\delta g} - \bar{F}_{\delta g} \rho]^T C_{\delta g \delta g}^{-1} [y_{\delta g} - \bar{F}_{\delta g} \rho] + [y_B - \bar{F}_B \chi]^T C_{BB}^{-1} [y_B - \bar{F}_B \chi] \quad (4) $$

where $\bar{F}_{\delta g}$ is a $m_{\delta g} \times n$ matrix, with $m_{\delta g}$ number of gravity observations and $n$ number of discretized densities, containing at the element $i, j$ the gravitational effect of the prism $j$ to the observation point $i$ (obtained by integrating $F_{\delta g}$). Similarly, $\bar{F}_B$ is a $m_B \times n$ matrix, obtained by integrating $F_B$. Eq. 4 represents the base of every joint inverse problem. At this point, we are left with a set of given observations $y_{\delta g}$, $y_B$ to be inverted in order to estimate a finite set of densities ρ and magnetic susceptibilities χ, within the discetized study region V. In such a problem we can consider as known the observations as well as their accuracy and spatial resolution. On the contrary, the volume V (defined by the coordinates of its vertices) and the size of the prisms in which it is discretized (namely Δx, Δy and Δz) should be defined before the inversion. In addition, to avoid edge effects, an external border region of side Δb should be also considered.

An important remark has to be made on the forward operators defined in Eq. 4. When dealing with such regional inversion problems, the gravitational effect of the Earth 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 and Capponi, 2019). Choosing the latter possibility, the first term of the data-misfit equation can be written as:

$$ \phi_{\delta g} = [y_{\delta g} - E[y_{\delta g}] - \left(I - \frac{J}{m_{\delta g}}\right) \bar{F}_{\delta g} \rho]^T C_{\delta g \delta g}^{-1} [y_{\delta g} - E[y_{\delta g}] - \left(I - \frac{J}{m}\right) \bar{F}_{\delta g} \rho] \quad (5) $$

where E[y<sub>δg</sub>] is the expectation operator applied to y<sub>δg</sub>, $m_{\delta g}$ is the number of observations, I is the identity matrix of size $m_{\delta g}$ and J is an $m_{\delta g} \times m_{\delta g}$ all-ones matrix. Basically $(I - \frac{J}{m_{\delta g}}) \bar{F}_{\delta g}$ represents a zero-mean forward operator (Sansó and Sampietro, 2022).

A similar argument can be used also for the magnetic field (at least for regions smaller than 2500 km × 2500 km). In fact, when trying to recover information on the Earth crust from the magnetic field, the lithospheric contribute 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., 2016, 2021). Total magnetic internal field models show a break in the slope of degree variances between degrees 13 and 16 (Langel and Estes, 1982) roughly corresponding to wavelengths between 3000 km 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 happen 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 implies 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 parametrization of the inverse problem, we will be lead 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 affect 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 fields 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 $(I - \frac{J}{m\delta_g})\bar{F}_{\delta_g}$ and $(I - \frac{J}{mB})\bar{F}_B$ operators. To evaluate these model errors, the effect of the studied volume (in terms of the grav/mag potential field) will be firstly 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.

## 3 | 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 firstly derive the accuracy of available gravity and magnetic models and then present the realistic synthetic model used to setup 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 who delivered the data, while for shipborne data available in GReD, having also raw observations, we inferred it from the analysis of crossovers following to 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 of second radial derivative (see e.g. Hofmann-Wellenhof and 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
<figure><img src="image_23.png" /><figcaption>FIGURE 1 XGM2019e error covariance in the Mediterranean Sea area</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 2 EMAG2v3, WDMAM-2 signal and error covariances.</figcaption></figure>

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

## 3.1 | Realistic Synthetic model

We then built a first synthetic model on the Mediterranean Sea area based on the following discontinuity surfaces:

1. bathymetry;

2. base of Plio-quaternary sediments;
<figure><img src="image_25.png" /><figcaption>FIGURE 3 Covariance functions of the reference geological horizons</figcaption></figure>

<table><caption>TABLE 1 Statistics of the reference geological horizons.</caption><thead><tr><th rowspan="2">LAYER</th><th>AVERAGE DEPTH</th><th>DEPTH STD</th><th>CORRELATION LENGTH</th></tr><tr><th>[m]</th><th>[m]</th><th>[m]</th></tr></thead><tbody><tr><td>Bathymetry</td><td>645</td><td>1090</td><td>15400</td></tr><tr><td>Base Plio-Quaternary</td><td>845</td><td>1274</td><td>18700</td></tr><tr><td>Base Mess. Salt</td><td>947</td><td>1647</td><td>17600</td></tr><tr><td>Basement</td><td>3943</td><td>3790</td><td>28600</td></tr><tr><td>Moho</td><td>31218</td><td>6810</td><td>33000</td></tr></tbody></table>

3. base of Messinian sediments;

4. basement;

5. Moho;

The bathymetry is taken from the GEBCO model (Mayer et al., 2018), the base of Plio-Quaternary and Messinian sediments have been derived by digitizing available models and seismic profiles (Capponi et al., 2020; Haq et al., 2020; de Voogd et al., 1992; Longacre et al., 2007; Feld et al., 2017; Makris and 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 NOAA (Straume et al., 2019), which however still largely depends in the Mediterranean Sea area from the CRUST2.0 global dataset (Molinari and Morelli, 2011). Finally we took the Moho depth from the Moho Map of the European Plate (Grad et al., 2009). In absence of better information we assume the Curie isotherm depth to be equal to the Moho discontinuity. Characteristics for the layers are reported in Tab. 1 and the respective covariance functions are reported in Fig. 3. The reference densities and magnetic susceptibilities are reported in Tab. 2. The crustal and upper mantle density models are taken from Fullea et al. (2021), while the susceptibility is derived from Hunt et al.
<table><caption>**TABLE 2** reference density and susceptibility model. * Crust and Upper mantle models lateral density and susceptibility distribution have been taken from Fullea et al. (2021) and Hunt et al. (1995).</caption><thead><tr><th>LAYER</th><th>AVERAGE DENSITY<br>[kg/m³]</th><th>DENSITY GRADIENT<br>[kg/m³/km]</th><th>SUSCEPTIBILITY<br>[10<sup>-6</sup>SI]</th></tr></thead><tbody><tr><td>Water</td><td>1030</td><td>0</td><td>-13</td></tr><tr><td>Plio-Quaternary</td><td>2220</td><td>0</td><td>10000</td></tr><tr><td>Messinian</td><td>2170</td><td>0</td><td>-30</td></tr><tr><td>Pre-Messinian</td><td>2400</td><td>6</td><td>1500</td></tr><tr><td>Crust</td><td>2890*</td><td>6</td><td>36796*</td></tr><tr><td>Upper Mantle</td><td>3321*</td><td>0</td><td>0</td></tr></tbody></table>

(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 3D model with spatial resolution of 1-arcminute (about 1.6 km) in x and y, 100 m in the z direction up to a depth of 50 km for a total of about 1.14 · 10⁹ cells. The expected accuracy of the shallower layers (base of Plio-quaternary and base of Messinian sediments) is of few hundreds metres while that of basement and moho is of few kilometres (as indicated in the literature used to build the model). Comparing these numbers with the STD of the amplitude of the different layers reported in Tab. 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 global level of these parameters for different sedimentary and crustal rocks (see e.g. Tenzer and Gladkikh (2014); Christensen and 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 vary also of order 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 previous published studies.

# 4 | RESULTS

## 4.1 | Selection of observations ($\delta g$, $T_{rr}$, dB)

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 ($\delta g$) and the second radial derivative of the potential ($T_{rr}$), both computed at sea level; for the magnetic field, the total magnetic intensity dB at 4 km above the surface. Fig. 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 Fig. 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 km 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
<figure><img src="image_26.png" /></figure>

**FIGURE 4** Computed gravitational effect of the synthetic model in terms of $\delta g$ (upper left panel) and $T_{rr}$ (lower left panel) and observed gravity field in terms of $\delta g$ (upper right panel) and $T_{rr}$ (lower right panel)

<figure><img src="image_27.png" /></figure>

**FIGURE 5** Computed total magnetic intensity anomaly due to the synthetic model (left panel), and observed total magnetic intensity (right panel).
<table><caption>TABLE 3 Signal-to-noise ratio for the different observables in the whole Mediterranean Sea area</caption><thead><tr><th>LAYER</th><th>δg</th><th>Trr</th><th>dF</th></tr></thead><tbody><tr><td>Water</td><td>29</td><td>4</td><td>2·10<sup>-4</sup></td></tr><tr><td>Plio-Quaternary</td><td>10</td><td>4</td><td>0.1</td></tr><tr><td>Messinian</td><td>14</td><td>2</td><td>3·10<sup>-4</sup></td></tr><tr><td>Pre-Messinian</td><td>90</td><td>8</td><td>0.04</td></tr><tr><td>Crust</td><td>252</td><td>11</td><td>2</td></tr><tr><td>Upper Mantle</td><td>203</td><td>5</td><td>0</td></tr></tbody></table>

that the computed forwards have 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 $T_{rr}$, $\delta g$ and $dF$. The signal-to-noise ratio, i.e. 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 second radial derivative (having an 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 $T_{rr}$ could be larger than the one of $\delta g$. Being interested in the whole Mediterranean Sea area our analysis will be performed assuming $\delta g$ as gravity observations. From Tab. 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 reaches 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.

## 4.2 | Planar vs. 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 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 spherical approximation of a given volume is performed by means of tesseract (Uieda et al., 2016; Baykiev 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 approximation for $\delta g$ and $T_{rr}$. Since, as stated in Sec. 2, 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 border effect, clearly visible with typical bell-shaped anomalies. These effects are
<figure><img src="image_28.png" /><figcaption>FIGURE 6 Difference between the gravitational effect of the reference volume in planar and spherical approximations in terms of $\delta g$ (left) and $T_{rr}$ (right).</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>FIGURE 7 Difference between the gravitational effect of the reference volume in planar and spherical approximations in terms of $\delta g$ (left) and $T_{rr}$ (right) once border effects are removed.</figcaption></figure>

easily removed both in planar and spherical approximation by assuming the borders to have a constant density $\bar{\rho}$:

* in planar approximation, this is done by computing the effect of a Bouguer slab with density $\bar{\rho}$ and thickness equal to 50 km and removing the effect of a prism with the same side of our reference volume and density $\bar{\rho}$.

* in spherical approximation the Bouguer slab is substituted by a spherical shell and the prism by a tesseract.

Once removed the border effect (here we used as $\bar{\rho}$ the average density of the synthetic model, i.e. 2938 kg/m³), we obtain the results shown in Fig. 7. We notice that the planar approximation in terms of $\delta g$ 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 case of planar approximation when using $\delta g$ 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 applied also to estimate the planar approximation error for the magnetic field. The difference between planar and spherical approximation, 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 Fig. 8. Similar results are obtained also for different directions of the inducing field as far as they are correctly projected between the spherical and planar reference systems.
<figure><img src="image_30.png" /><figcaption>**FIGURE 8** Difference between the Total magnetic intensity of the reference volume in planar and spherical approximations.</figcaption></figure>

## 4.3 | 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 Fig. 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, i.e. 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 and 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 a 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 tesseracts for the spherical approximation and prism for the planar one. We show in Fig. 9 the results obtained for the planar approximation in terms of δg, $T_{rr}$ and dF. Starting from the gravity field, we observe that a border of about 15° should be consider in order to reduce edge effects STD to values lower than the observation error. As expected this value is notably reduced when moving to $T_{rr}$ to less than 5°. A behaviour similar to the one of $T_{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 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 tesseracts to compute the radial component of the crustal magnetic field on two areas for different border size. The situation is very similar also when moving to the spherical approximation, with δg 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 δg and 10° for $T_{rr}$ is required.
<figure><img src="image_31.png" /><figcaption>**FIGURE 9** Error (border effect) as a function of the size of the border for $\delta g$ (upper panel), $T_{rr}$ (middle panel) and $dF$ (lower panel).</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>FIGURE 10 Error (mantle effect) as a function of the maximum depth in terms of $\delta g$.</figcaption></figure>

## 4.4 | Maximum depth of the 3D model

Similarly to the border size in the latitude, longitude (North-South, West-East) directions, in order to deal with gravity we also have to analyse the maximum required depth of the considered 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 contribute) is expected. Supposing that the mantle density distribution is known with sufficient accuracy (see Sampietro and 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 details 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 *z* direction of 1 km. In Fig. 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 km 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 399 km and 400 km has an STD of about 0.01 mGal, thus making the effect of the deeper density distribution negligible.

## 4.5 | Model resolution

We use our synthetic model also 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 *n*. We increase *n* as long as the downsampling error is smaller than the observation one. We report in Fig. 11 maps of the downsampling error in terms of $\delta g$ for *n* = 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). Fig. 12 shows the STD of this downsampling error in the *x* and *y* directions. In the figure the observation error for $\delta g$ and $T_{rr}$ is also reported. We see how, considering the whole Mediterranean Sea area and relying on the XGM2019e and the EMAG2v3 models *n* = 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 *n* = 20 will reduce the total number of cells of our 3D model from about $1.14 \cdot 10^9$ to about $2.84 \cdot 10^6$ cells. When we
<figure><img src="image_33.png" /><figcaption>FIGURE 11 Downsampling error in terms of $\delta g$ for $n = 2, 5, 20, 40$</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>FIGURE 12 Error from downsampling in the x and y directions for the considered functionals. Dashed lines represent Observation errors for $\delta g$ (blue) and $T_{rr}$. Units are mGal, E, nT for $\delta g$, $T_{rr}$ and dF respectively</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>FIGURE 13 Error due to the downsampling in the *z* direction for the considered functionals. Dashed lines represent observation errors for *δg* and *Trr*. Units are mGal, E, nT for *δg*, *Trr* and dF respectively</figcaption></figure>

move to higher values of *n* the total error start to be dominated by the model one thus 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 performed also 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 model. Results are reported in Fig. 13 where it can be seen that a resolution of about 600 m is enough to discretize the volume in the vertical direction.

5 | 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, in 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 non-uniqueness 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 simplifying the volume to a two-layer problem, i.e. 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 and 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³ for the crust and of 3330 kg/m³ for the mantle separated by the Moho model presented in the Sec. 3. 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 been also simulated. The gravitational effect due to the true model and the simulated observation error are reported in Fig. 14. As 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
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

<figure><img src="image_36.png" /><figcaption>FIGURE 14 gravitational effect due to the mass distribution of the true model (left panel) and simulated observation error (right panel).</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>FIGURE 15 True model (up), and a-priori information for the inversion (down). Green lines represent the range of variability of the a-priori Moho.</figcaption></figure>

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 Fig. 15.

We start by evaluating the effect of the border size. This is done by performing the inversion with borders of 4000 km, 1500 km, 1000 km 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 Fig 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 Fig. 16 are basically the combined effect of observation errors and of errors due to the specific inversion used. Fig. 17 shows the differences between the true and estimated Moho for the other border sizes. For border 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 wavelengths 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
<figure><img src="image_38.png" /><figcaption>FIGURE 16 Difference between the true and estimated Moho with 4000 km border.</figcaption></figure>

improve the inversion result.

In a second test we consider the model spatial resolution in the planar direction. We perform different inversions with spatial resolutions of 8 km, 33 km and 66 km. Looking to the differences with respect to the true model (Fig. 18), we see that the 8 km and 33 km spatial resolutions show almost the same error with STD of 394 m and 402. 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 concentrated in the Black Sea, close to Sicily or between Sardinia and Balearic Islands. We notice that in terms of cells number we move from about $8 \cdot 10^6$ cells of the 8 km spatial resolution model to the $5 \cdot 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.

## 6 | 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, computed with a set of forwards performed on a realistic synthetic model of the considered volume.

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 future study dealing with regional areas. Looking to the results in the Mediterranean Sea region we firstly assessed the accuracy of global gravity and magnetic field models: the XGM2019e global gravity field model in terms of $\delta g$ synthesised 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 3D model in terms of geological units, densities and susceptibilities on the considered volume. Relying on the above observations errors and synthetic model we performed a set of simulation to prove that even for such a large region the planar approximation is sufficient entailing errors smaller than the observations ones. Looking to the signal-to-noise ratio for the different geological horizons of the synthetic model we also observed that, on the
<figure><img src="image_39.png" /><figcaption>**FIGURE 17** Differences between the true and reference Moho for different border sizes.</figcaption></figure>
<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /><figcaption>FIGURE 18 Differences between the true and reference Moho for different spatial resolutions.</figcaption></figure>
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 x and y 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 tesseract. 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 hours is required. However the application of well known acceleration by means of 1D FFT and the use of parallel computing are expected to reduce the computational time of 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 right rectangular prism and FFT requires less than 5 seconds, thus allowing to set all the parameters investigated in less than 10 minutes.

Some inversions tests, performed on a synthetic case study, gives 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 a proper knowledge of the volume investigated and of the error of the available observations. The former is usually available at least in a rough way from 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 literature review and knowing how the observations have been acquired and processed.

## Acknowledgements

Acknowledgements should include contributions from anyone who does not meet the criteria for authorship (for example, to recognize contributions from people who provided technical help, collation of data, writing assistance, acquisition of funding, or a department chairperson who provided general support), as well as any funding or other support information.

## Conflict of Interest

Daniele Sampietro, Martina Capponi, Erwan Thebauld, Lydie Gailler declare that they have no conflict of interest.

## Data Availability Statement

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
references

An, M. (2012) A simple method for determining the spatial resolution of a general inverse problem. *Geophysical Journal International*, **191**, 849–864.

Backus, G., George, B., Parker, R. L., Parker, R. and Constable, C. (1996) *Foundations of geomagnetism*. Cambridge University Press.

Baykiev, E., Ebbing, J., Brönner, M. and Fabian, K. (2016) Forward modeling magnetic fields of induced and remanent magnetization in the lithosphere using tesseract. *Computers & geosciences*, **96**, 124–135.

Bhattacharyya, B. (1964) Magnetic anomalies due to prism-shaped bodies with arbitrary polarization. *Geophysics*, **29**, 517–531.

Blakely, R. J. (1996) *Potential theory in gravity and magnetic applications*. Cambridge university press.

Bosch, M., Meza, R., Jiménez, R. and Hönig, A. (2006) Joint gravity and magnetic inversion in 3d using monte carlo methods. *Geophysics*, **71**, G153–G156.

Capponi, M., GlaviCh, E. and Sampietro, D. (2020) Map of plio-quaternary sediment depths in the mediterranean sea. *Bollettino di Geofisica Teorica e Applicata*, **61**, 421–432.

Christensen, N. I. and Mooney, W. D. (1995) Seismic velocity structure and composition of the continental crust: A global view. *Journal of Geophysical Research: Solid Earth*, **100**, 9761–9788.

Feld, C., Mechie, J., Hübscher, C., Hall, J., Nicolaides, S., Gurbuz, C., Bauer, K., Louden, K. and Weber, M. (2017) Crustal structure of the eratosthenes seamount, cyprus and s. turkey from an amphibian wide-angle seismic profile. *Tectonophysics*, **700**, 32–59.

Fullea, J., Lebedev, S., Martinec, Z. and Celli, N. (2021) Winterc-g: mapping the upper mantle thermochemical heterogeneity from coupled geophysical-petrological inversion of seismic waveforms, heat flow, surface elevation and gravity satellite data. *Geophysical Journal International*, **226**, 146–191.

Galley, C. G., Lelièvre, P. G. and Farquharson, C. G. (2020) Geophysical inversion for 3d contact surface geometry. *Geophysics*, **85**, K27–K45.

Grad, M., Tiira, T. and Group, E. W. (2009) The moho depth map of the european plate. *Geophysical Journal International*, **176**, 279–292.

Haq, B., Gorini, C., Baur, J., Moneron, J. and Rubino, J.-L. (2020) Deep mediterranean's messinian evaporite giant: How much salt? *Global and Planetary Change*, **184**, 103052.

Hofmann-Wellenhof, B. and Moritz, H. (2006) *Physical geodesy*. Springer Science & Business Media.

Hunt, C. P., Moskowitz, B. M., Banerjee, S. K. et al. (1995) Magnetic properties of rocks and minerals. *Rock physics and phase relations: A handbook of physical constants*, **3**, 189–204.

Kuhn, M., Featherstone, W. and Kirby, J. (2009) Complete spherical bouguer gravity anomalies over australia. *Australian journal of Earth sciences*, **56**, 213–223.

Langel, R. and Estes, R. (1982) A geomagnetic field spectrum. *Geophysical Research Letters*, **9**, 250–253.

Laske, G., Masters, G., Ma, Z. and Pasyanos, M. (2013) Update on crust1.0—a 1-degree global model of earth's crust. In *Geophys. Res. Abstr*, vol. 15, 2658.

Lesur, V., Hamoudi, M., Choi, Y., Dyment, J. and Thébault, E. (2016) Building the second version of the world digital magnetic anomaly map (wdmam). *Earth, Planets and Space*, **68**, 1–13.
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

Li, X. (2001) Vertical resolution: Gravity versus vertical gravity gradient. *The Leading Edge*, **20**, 901–904.

Li, Y. and Oldenburg, D. W. (1996) 3-d inversion of magnetic data. *Geophysics*, **61**, 394–408.

– (1998a) 3-d inversion of gravity data. *Geophysics*, **63**, 109–119.

– (1998b) Separation of regional and residual magnetic field data. *Geophysics*, **63**, 431–439.

Longacre, M., Bentham, P., Hanbal, I., Cotton, J. and Edwards, R. (2007) New crustal structure of the eastern mediterranean basin: detailed integration and modeling of gravity, magnetic, seismic refraction, and seismic reflection data. In *EGM 2007 International Workshop: Innovation in EM, Grav and Mag Methods: A New Perspective for Exploration*, vol. 15, 18.

Makris, J. and Yegorova, T. (2006) A 3-d density–velocity model between the cretan sea and libya. *Tectonophysics*, **417**, 201–220.

Mayer, L., Jakobsson, M., Allen, G., Dorschel, B., Falconer, R., Ferrini, V., Lamarche, G., Snaith, H. and Weatherall, P. (2018) The nippon foundation–gebco seabed 2030 project: The quest to see the world’s oceans completely mapped by 2030. *Geosciences*, **8**, 63.

Meyer, B., Saltus, R. W. and Chulliat, A. (2016) Emag2-v3: A new global compilation of lithospheric magnetic anomalies. In *AGU Fall Meeting Abstracts*, vol. 2016, GP43A–1210.

Molinari, I. and Morelli, A. (2011) Epcrust: a reference crustal model for the european plate. *Geophysical Journal International*, **185**, 352–364.

Moritz, H. (1990) The figure of the earth: theoretical geodesy and the earth’s interior. *Karlsruhe: Wichmann*.

Murböck, M. (2015) *Virtual constellations of next generation gravity missions*. Ph.D. thesis, Technische Universität München.

Nabighian, M. N., Grauch, V., Hansen, R., LaFehr, T., Li, Y., Peirce, J. W., Phillips, J. D. and Ruder, M. (2005a) The historical development of the magnetic method in exploration. *Geophysics*, **70**, 33ND–61ND.

– (2005b) The historical development of the magnetic method in exploration. *Geophysics*, **70**, 33ND–61ND.

Nagy, D. (1966) The gravitational attraction of a right rectangular prism. *Geophysics*, **31**, 362–371.

Oldenburg, D. W. (1974) The inversion and interpretation of gravity anomalies. *Geophysics*, **39**, 526–536.

Pasyanos, M. E., Masters, T. G., Laske, G. and Ma, Z. (2014) Litho1. 0: An updated crust and lithospheric model of the earth. *Journal of Geophysical Research: Solid Earth*, **119**, 2153–2173.

Sampietro, D. (2015) Geological units and moho depth determination in the western balkans exploiting goce data. *Geophysical Journal International*, **202**, 1054–1063.

Sampietro, D. and Capponi, M. (2019) Practical tips for 3d regional gravity inversion. *Geosciences*, **9**, 351.

Sampietro, D., Capponi, M., Mansi, A., Gatti, A., Marchetti, P. and Sansò, F. (2017) Space-wise approach for airborne gravity data modelling. *Journal of Geodesy*, **91**, 535–545.

Sampietro, D. and Sansò, F. (2012) Uniqueness theorems for inverse gravimetric problems. In *VII Hotine-Marussi Symposium on Mathematical Geodesy*, 111–115. Springer.

Sansó, F. and Sampietro, D. (2022) Analysis of the gravity field.

Straume, E. O., Gaina, C., Medvedev, S., Hochmuth, K., Gohl, K., Whittaker, J. M., Abdul Fattah, R., Doornenbal, J. C. and Hopper, J. R. (2019) Globsed: Updated total sediment thickness in the world’s oceans. *Geochemistry, Geophysics, Geosystems*, **20**, 1756–1772.
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

Tenzer, R. and Gladkikh, V. (2014) Assessment of density variations of marine sediments with ocean and sediment depths. *The Scientific World Journal*, **2014**.

Thébault, E., Hulot, G., Langlais, B. and Vigneron, P. (2021) A spherical harmonic model of earth's lithospheric magnetic field up to degree 1050. *Geophysical Research Letters*, **48**, e2021GL095147.

Thébault, E., Vigneron, P., Langlais, B. and Hulot, G. (2016) A swarm lithospheric magnetic field model to sh degree 80. *Earth, Planets and Space*, **68**, 1–13.

Uieda, L., Barbosa, V. C. and Braitenberg (2016) Tesseracts: Forward-modeling gravitational fields in spherical coordinates. *Geophysics*, **81**, F41–F48.

de Voogd, B., Truffert, C., Chamot-Rooke, N., Huchon, P., Lallemant, S. and Le Pichon, X. (1992) Two-ship deep seismic soundings in the basins of the eastern mediterranean sea (pasiphae cruise). *Geophysical Journal International*, **109**, 536–552.

Zahorec, P., Papčo, J., Pašteka, R., Bielik, M., Bonvalot, S., Braitenberg, C., Ebbing, J., Gabriel, G., Gosar, A., Grand, A. et al. (2021) The first pan-alpine surface-gravity database, a modern compilation that crosses frontiers. *Earth System Science Data*, **13**, 2165–2209.

Zeng, X., Wan, X., Lin, M. and Wang, W. (2022) Gravity field forward modelling using tesseracts accelerated by taylor series expansion and symmetry relations. *Geophysical Journal International*, **230**, 1565–1584.

Zingerle, P., Pail, R., Gruber, T. and Oikonomidou, X. (2020) The combined global gravity field model xgm2019e. *Journal of Geodesy*, **94**, 1–12.