# From space to lithosphere: inversion of the GOCE gravity gradients. Supply to the Earth's interior study

Matthieu Plasman, Christel Tiberi, Cecilia Cadio, Anita Thea Saraswati, Gwendoline Pajot-Métivier, Michel Diament

► To cite this version:

Matthieu Plasman, Christel Tiberi, Cecilia Cadio, Anita Thea Saraswati, Gwendoline Pajot-Métivier, et al.. From space to lithosphere: inversion of the GOCE gravity gradients. Supply to the Earth's interior study. Geophysical Journal International, 2020, 223 (1), pp.398-419. 10.1093/gji/ggaa318. hal-03010550

HAL Id: hal-03010550

https://hal.science/hal-03010550

Submitted on 11 Jun 2021

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# From space to lithosphere: inversion of the GOCE gravity gradients. Supply to the Earth's interior study

Matthieu Plasman, Christel Tiberi, Cecilia Cadio, Anita Saraswati, Gwendoline Pajot-Métivier, Michel Diament

► To cite this version:

Matthieu Plasman, Christel Tiberi, Cecilia Cadio, Anita Saraswati, Gwendoline Pajot-Métivier, et al. From space to lithosphere: inversion of the GOCE gravity gradients. Supply to the Earth's interior study. Geophysical Journal International, Oxford University Press (OUP), 2020, 10.1093/gji/ggaa318. hal-03012530

HAL Id: hal-03012530

https://hal.archives-ouvertes.fr/hal-03012530

Submitted on 18 Nov 2020

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.
**From space to lithosphere : inversion of the GOCE gravity gradients. Supply to the Earth's interior study.**

Matthieu Plasman<fnref n="1" />,<fnref n="3" />, Christel Tiberi<fnref n="1" />, Cecilia Cadio<fnref n="1" />, Anita Thea Saraswati<fnref n="2" />,
Gwendoline Pajot-Métivier<fnref n="3" />,<fnref n="4" /> and Michel Diament<fnref n="3" />

**SUMMARY**

The emergence of high resolution satellite measurements of the gravitational field (GOCE mission) offers promising perspectives for the study of the Earth's interior. These new data call for the development of innovative analysis and interpretation methods. Here we combine a forward prism computation with a Bayesian resolution approach to invert for these gravity gradient data configuration. We apply and test our new method on satellite data configuration, i.e. 225 km height with a global and homogeneous geographic distribution. We first quantify the resolution of our method according to both data and parameterization characteristics. It appears that for reasonable density contrast values (0.1 g.cm⁻<fnref n="3" />) crustal structures have to be wider than ~28 km to be detectable in the GOCE signal. Deeper bodies are distinguishable for greater size (35 km size at 50 km depth, ~80 km at 300 km depth). We invert the six tensor components, among which five are independent. By carefully testing each of them and their different combinations, we enlighten a trade off between the recovery of data and the sensitivity to inversion parameters. We particularly discussed this characteristic in terms of geometry of the synthetic model tested

<footnote n="1">¹ Géosciences Montpellier, CNRS-UMR5243, Université de Montpellier, 34000 Montpellier, France</footnote>

<footnote n="2">² Faculté des Sciences, de la Technologie et de la Communication, Université du Luxembourg, Esch-sur-Alzette, Luxembourg</footnote>

<footnote n="3">³ Université de Paris, Institut de physique du globe de Paris, CNRS, IGN, F-75005 Paris, France</footnote>

<footnote n="4">⁴ ENSG-Géomatique, IGN, F-77455 Marne-la-Valle, France</footnote>
(structures orientation, 3-D geometry, etc.). In terms of RMS value, each component is
always better explained if inverted solely, but the result is strongly affected by the in-
version parameterization (smoothing, variances, etc.). On the contrary, the simultaneous
inversion of several components displays a significant improvement for the global ten-
sor recovery, more dependent on data than on density variance or on smoothness control.
Comparing gravity and gradient inversions, we highlight the superiority of the GG data to
better reproduce the structures especially in terms of vertical location. We successfully test
our method on a realistic case of a complex subduction case for both gradient and gravity
data. While the imaging of small crustal structures requires terrestrial gravity dataset, the
longest wavelength of the slab is well recovered with both data sets. The precision and
homogeneous coverage of GOCE data however, counterbalance the heterogeneous and
often quite nonexistence coverage of terrestrial gravity data. This is particularly true in
large areas which requires a coherent assemblage of heterogeneous data sets, or in high
relief, vegetally covered and offshore zones.

**Key words:** GOCE satellite, Full Gradient Tensor, Gravity Gradients Data, Gravity method, Geophysical inversion

# 1 INTRODUCTION

The recent satellite mission GOCE (Gravity field and steady-state Ocean Circulation Explorer) has opened up the scope for new kind of gravitational field observations by measuring globally for the first time the full gravity gradient tensor. With its three pairs of accelerometers sensitive to the gravitational acceleration in the three directions of the satellite reference system, GOCE measures the second order derivatives of the gravitational potential in three independent directions (Drinkwater et al. 2006). Initially deployed on an orbit at 255 km height, it was then reduced to 225 km in its lower orbit phase at the end of the mission to reach higher precision and resolution. In this way, GOCE provides homogeneous global coverage measurements with a resolution of 80 km for a precision better than 1 mGal/10 mE at the sea level (Rummel et al. 2011; Bouman et al. 2015). The main objectives of the space mission were the study of the mean sea level (Knudsen et al. 2011), the thermohaline circulation (Janjić et al. 2012) and the ice thickness (Hirt 2014) to quantify the climate change and its impacts through the measurement of the Earth gravity field. Moreover, the accurate measurement of the gravity field and the geoid at this scale bring some new and precious geophysical informations on e.g. the
Earth geodynamics (Braitenberg et al. 2011; Panet et al. 2014; Braitenberg 2015), mantle composition (Cadio et al. 2011) or the geometry of the lithosphere (Reguzzoni & Sampietro 2012; van der Meijde et al. 2015a).

In the study of the Earth's interior, the understanding of the interactions between local and global processes is a fundamental step to apprehend its structure and global dynamics. The knowledge of the complete gravitational field is then essential to tackle this question. Ground gravity measurements are since many decades and the gravity method is now a classical approach. As a consequence and since several years, the airborne approach, with gravity or gravity gradient measurements, has become more widely developed for geophysical exploration (Dransfield & Milkereit 2007), especially in the industrial prospecting (mining companies). It has led to the creation of new and more precise instruments (e.g. GREMLIT project, Douch et al. 2014), analysis methods (e.g. Forsberg & Olesen 2010; Hassan 2005) and inversion schemes (e.g. Martinez et al. 2012; Geng et al. 2014). The development of new methods based on the inversion of the whole gravity gradient tensor at satellite height is now the consequent step to fully complete the informations brought by the ground and airborne gravity gradient data.

Compared to gravity methods which only use the vertical component of $\vec{g}$ vector, the gravity gradient (GG) method manages 5 independent components of the gravity gradient tensor ($T_{xx}, T_{xy}, T_{xz}, T_{yy}$ and $T_{yz}$). It results in a better resolving capability and a higher frequency content (Beiki 2010). The inversion of gravity data ($g_z$) is classically used in geophysics either with terrestrial (e.g. Jacoby & Smilde 2009), airborne or satellite (e.g. Cambiotti & Sabadini 2012) measurements. Even if the GG method is more and more used, particularly in the oil or mineral domain and/or for crustal targets (Martinez et al. 2012), only few studies tried to invert the GG GOCE data so far, particularly in lithospheric studies (e.g. Afonso et al. 2019). First attempts to invert GOCE data were based on the gravity anomalies deduced from the GG data and did not consider the full tensor gradient (FTG) components (Basuyau et al. 2013; O'Donnell et al. 2013; Dogru et al. 2018). Only recently a few studies have directly used the GG data from GOCE in an inversion process (e.g. Li et al. 2017).

However the benefit of gravity gradient measurements in an inversion process, compared to the gravity data, can be reasonably questioned. Theoretically the derivatives of the gravity field do not contain more information than the gravity field itself (Pedersen & Rasmussen 1990). With a perfect knowledge of the potential at each point we would also know its derivatives without measuring them. But in practice the precision and heterogeneous data coverage prevents from perfectly knowing it in each point. Moreover it will be necessary to make a data prolongation to compute the complete tensor (Pedersen & Rasmussen 1990). The redundancy (six components) and complementarity (each component brings information in a specific direction) of the measured GG data allow to reduce the data
noise and to provide a more complete directional information to better locate the sources (Pajot 2007).
However for near surface data, Bell et al. (1997) have shown that the joint interpretation of gravity and
gravity gradient data is needed to improve the resolution of the small wavelength structures. It is why
the gravity gradient measurement and its inversion appear essential in geophysical imaging.

In this study we present a new methodology for GOCE gravity gradient data inversion based on
the discretization of the Earth's interior into 3-D rectangular prism. The gravity gradient effect of
each prism is computed with the Uieda et al. (2016) method. The inversion can be resumed to the
minimization of the difference between the observed and the synthetic data to retrieved the best 3-D
distribution of rectangular prisms. This minimization is then realized with a Bayesian optimization
(Menke 1984) to include a priori and a posteriori informations. An originality of our approach is to
considere in a joint scheme all the gradient components, or a combination of them, or a single one.

We simultaneously invert the full tensor components ($T_{xx}$, $T_{xy}$, $T_{xz}$, $T_{yy}$, $T_{yz}$ and $T_{zz}$) to obtain
density contrast variations within the inquired volume. It offers the opportunity 1) to better constrain
the geometry of the structures and 2) to reduce the non-uniqueness of the inversion, which is an impor-
tant limitation of the gravity inversion. Moreover, compared to existing terrestrial gravity data, GOCE
data display a better sensitivity to deeper and larger structures thanks to global and homogeneous data
coverage. It thus complements other sources of gravity or gradient data. However the contribution of
GG data in an inversion process needs to be quantified compared to the gravity inversion, first with
the GOCE data configuration and then in a more general case. In order to fully evaluate the resolu-
tion power of our method, we focus in this paper on synthetic tests. It allows the exemption of data
uncertainty and the combined effect of several sources.

## 2 GOCE DATA IN THE EARTH'S INTERIOR STUDY: STATE OF THE ART

The use of the GOCE gravity gradient in an inversion scheme requires some considerations. The first
concern is of technical order. Indeed, because of the gradiometers configuration, the $T_{xx}$, $T_{yy}$, $T_{zz}$ and
$T_{xz}$ are determined with high accuracy whereas $T_{xy}$ and $T_{yz}$ are around twice less accurate (Bouman
et al. 2015). Moreover the data are first available in the GRF (Gradiometer Reference Frame) which
co-rotates with the satellite (Bouman et al. 2015) making their use quite complex. Classically, people
prefer to use the Local North-Oriented Frame (LNOF) (Bouman et al. 2009, 2011; Fuchs & Bouman
2011), a kind of compromise between the original GOCE data (given along the orbit) and an easy use.
Finally an alternative representation exists with the GG grid at the mean satellite altitude (Bouman
et al. 2015).

The second consideration is, independently of the frame of the data, that the gravity gradient maps
of the 6 components are very hard to interprete as soon as we are far from a simple synthetic case (e.g. a
prism). Several methods have been developed to process and extract the maximal information of these
complex signals without any inversion procedure. For example Analytic signal technique, Euler or
Werner deconvolution (e.g. Beiki 2010; Beiki & Pedersen 2011; Bell et al. 1997; Li 2001; Mikhailov
et al. 2007; Pajot et al. 2008) help determining first order geographic location (horizontal position or
depth) and geometry (width, dip, etc.). However these methods are strongly dependent on a certain
number of hypothese and (over)simplifications (geometry approximation, ponctual source, constant
density, etc.). That is why they appear to be more preliminary methods that bring some constraints to
the inversion.

Very few inversion methods have been applied or developed for the GOCE GG, even though this
approach allows to pass through the apparent complexity of the signal to focus only on their meaning
by recovering the geometry and location of the sources. Two types of approaches have been applied.
One based on the inversion of the Bouguer anomaly, deduced from the GOCE data, in a "classic"
gravity inversion method (e.g. Basuyau et al. 2013). The other one, less widespread, directly inverts
the GG data. For example Li et al. (2017) use the six tensor components for their inversion, based on
the Tikhonov regularization theory. Recently, Tian et al. (2019) developed a preconditioned conjugate
gradient inversion approach for the GG data and applied it to synthetic (e.g. prisms) and real (North
China Craton) cases. In these studies, data are downwarded to 10 km above the ground surface. Even
if their results underlines the feasibility of GG inversion, this highlights 1) the difficulty to assess the
real benefit of the GOCE data and 2) the downward continuation may improve the total signal power
but the high frequency content increases the data noise as well as the near-surface density induced
signal (e.g. topography) which is not desirable for lithospheric studies (Bouman et al. 2015).

Several models of gravity field have been created from GOCE data. They can be GOCE only
model (Pail et al. 2011), combined models from satellite data only (Bruinsma 2010) or a mix of ground
and space data (Migliaccio et al. 2010; Pail et al. 2011). These combinations are aimed to produce the
most accurate gravity models. In a GOCE/GRACE comparison, Pail et al. (2011) showed that GOCE
displays better performance than GRACE (from degree 150 and beyond). These different data sets
have been used in different geodynamic contexts, to explore different structures and processes. We
can however retain, from shallow to deeper objects, the van der Meijde et al. (2015b) and Sampietro
et al. (2014) studies to image the Moho. Dogru et al. (2018) used the GOCE based gravity field model
to focus on the lithospheric structure (Moho, LAB depth) in the western Anatolia. Deeper, Panet et al.
(2014) imaged the mass distribution of the Earth's mantle.

A major interrogation remains about the choice of the data to be inverted. Is it better to use the
original measured GOCE gradient data (in the GRF or LNOF system) when the use of GG grids
seems to be more adapted to the inversion (homogeneous coverage, reduced data number, constant
height) ? It appears that these grids present three significant advantages. First, they prevent error from rotating frames at different altitudes and varying error characteristics (Bouman et al. 2016). Moreover, compared to global models, GG grids also show higher signal content due to a less regularized solution and are also closer to the observations measured at orbit height (Bouman et al. 2016). This is why we choose to use gridded synthetic gravity gradient data at the satellite height in our study.

## 3 INVERSION PROCEDURE

Our inversion method uses GG data to image the Earth's interior in terms of 3-D density patterns. The inspected volume is discretized into rectangular prisms. The inversion consists of a forward calculation to compute the GG response of all prisms, associated with a minimization algorithm that searches for the lowest misfit between observed and predicted GG data.

### 3.1 Forward modelling

#### 3.1.1 Data definition

The gravity gradient tensor contains the second derivatives of the Earth's gravitational potential (V) in the (x, y, z) directions. It corresponds to a second order tensor which, in a Cartesian coordinate system, can be written as :

$$ T = \begin{pmatrix} T_{xx} & T_{xy} & T_{xz} \\ T_{yx} & T_{yy} & T_{yz} \\ T_{zx} & T_{zy} & T_{zz} \end{pmatrix} = \begin{pmatrix} \frac{\partial^2 V}{\partial x^2} & \frac{\partial^2 V}{\partial y \partial x} & \frac{\partial^2 V}{\partial z \partial x} \\ \frac{\partial^2 V}{\partial x \partial y} & \frac{\partial^2 V}{\partial y^2} & \frac{\partial^2 V}{\partial z \partial y} \\ \frac{\partial^2 V}{\partial x \partial z} & \frac{\partial^2 V}{\partial y \partial z} & \frac{\partial^2 V}{\partial z^2} \end{pmatrix} \quad (1) $$

In any orthogonal coordinates system, T is a symmetric tensor ($T_{ij} = T_{ji}$). From Laplace equation the trace of the tensor is null, which leads to dependent $T_{xx}$, $T_{yy}$ and $T_{zz}$. Thus, when we use the six components $T_{xx}$, $T_{xy}$, $T_{xz}$, $T_{yy}$, $T_{yz}$ and $T_{zz}$ of the tensor, only 5 are mathematically independent. However because off-diagonal components are sensitive to variations of the gravitational potential in two different directions, the six tensor components display a kind of interdependency when characterizing a same object.

#### 3.1.2 Model parameterization and direct computation

The investigated volume is discretized in layers, subdivided into parallelepiped blocks, each associated with a constant density contrast value. The size of the prisms can vary within a layer, allowing to model irregular structures. The computation of the GG response for the whole volume is performed
by adding, for each tensor component, the gradient effect of each rectangular prism in every layers (Uieda et al. 2016).

Considering a 3-D parallelepiped (Fig. 1), with $x_1, y_1, z_1$ and $x_2, y_2, z_2$ the coordinates of its opposite corners, and according to the law of gravitation, the resulting gravity effect at a point $(x_0, y_0, z_0)$, outside this structure, is described by :

$$g_z(x_0, y_0, z_0) = G\rho \int_{x_1}^{x_2} \int_{y_1}^{y_2} \int_{z_1}^{z_2} \frac{z-z_0}{((x-x_0)+(y-y_0)+(z-z_0))^3} dxdydz \quad (2)$$

with G the gravitational constant and ρ the density.

An approximate expression of the gravity gradient components can be deduced from equation 2 (Nagy et al. 2000) :

$$T_{xx} = |||G \arctan \frac{(x - x_0)r}{(y - y_0)(z - z_0)}| \frac{x_2}{x_1} | \frac{y_2}{y_1} | \frac{z_2}{z_1} \rho \quad (3)$$

$$T_{yy} = |||G \arctan \frac{(y - y_0)r}{(x - x_0)(z - z_0)}| \frac{x_2}{x_1} | \frac{y_2}{y_1} | \frac{z_2}{z_1} \rho \quad (4)$$

$$T_{zz} = |||G \arctan \frac{(z - z_0)r}{(x - x_0)(y - y_0)}| \frac{x_2}{x_1} | \frac{y_2}{y_1} | \frac{z_2}{z_1} \rho \quad (5)$$

$$T_{xy} = |||G \ln(z - z_0 + r)| \frac{x_2}{x_1} | \frac{y_2}{y_1} | \frac{z_2}{z_1} \rho \quad (6)$$

$$T_{xz} = |||G \ln(y - y_0 + r)| \frac{x_2}{x_1} | \frac{y_2}{y_1} | \frac{z_2}{z_1} \rho \quad (7)$$

$$T_{yz} = |||G \ln(x - x_0 + r)| \frac{x_2}{x_1} | \frac{y_2}{y_1} | \frac{z_2}{z_1} \rho \quad (8)$$

With this approach (eqs 3-8), the second derivatives of $g_z$ (eq. 2) are not defined for points located inside or at the surface of a prism. But in our case, this configuration cannot arise, simply due to the data acquisition height (225 km) compared with the Earth surface.

We compute the effect of the rectangular prism using TESSEROID algorithm (Uieda et al. 2016).

## 3.2 Minimization approach

The inversion process consists of minimizing the difference between the calculated ($d_{cal}$) and the observed ($d_{obs}$) GG data. For our dedicated problem we use a $L_2$ norm (least squares). Geophysical inversion being a typically ill-posed problem, we use a Bayesian approach, that allows us to include a
priori and a posteriori information through the initial model and covariance matrices (Menke 1984). We suppose a Gaussian probability density function for both data and parameters distribution, and use the covariance matrices $C_d$ and $C_p$ to weight the data and the parameters distribution, respectively. We have then to minimize the following expressions for the data part:

$$ (d_{obs} - d_{cal})^T C_d^{-1} (d_{obs} - d_{cal}) \qquad (9) $$

and for the parameter :

$$ (p - p_0)^T C_p^{-1} (p - p_0) \qquad (10) $$

The symbol $^T$ denotes the transpose operation, $p_0$ is the density starting model and $p^T = (\Delta\rho)$ is the optimum density distribution. We include an additional smoothing constraint in our scheme to limit unrealistic spatial density changes, that could appear in gravity inversion:

$$ \left( \frac{\Delta p}{\Delta R} \right)^T C_s^{-1} \left( \frac{\Delta p}{\Delta R} \right) \qquad (11) $$

Equation 11 relates adjacent density contrasts $\Delta p$ with their corresponding distance $\Delta R$ (first derivatives). The matrix $C_s$ is diagonal and contains the smoothing factor, allowing to weight this supplementary condition.

Minimizing the sum of these three equations means to solve the following expression:

$$ p = p_0 + (A^t C_d^{-1} A + C_p^{-1} - C_s^{-1} D_s)^{-1} \cdot (A^t C_d^{-1} (d_{obs} - d_{calc}) + C_s^{-1} s) \qquad (12) $$

where $A$ is the partial derivatives matrix of the calculated data $d_{calc}$ relative to the parameters $p$. $D_s$ and $s$ are the matrix and the vector, respectively, that correspond to the roughness control of the model (smoothing). Their structure is as follows (for all blocks $j$ adjacent to block $i$ and $j \neq i$):

$$ D_{ii} = - \sum_{j} R_{ij}^{-2} \qquad (13) $$

$$ D_{ij} = R_{ij}^{-2} \qquad (14) $$

$$ s_i = \sum_j \frac{p_i - p_j}{R_{ij}^2} \qquad (15) $$

After linearization of eq. 12, we obtain an iterative procedure for the estimation of the optimum density model $p$. A new density model $p$ is computed with respect to the previous density distribution $p_0$. The convergence of the inversion is verified through the root mean square (RMS) decrease. The process is stopped either after a determined number of iterations, or when a given threshold ($\epsilon$) is
reached ((dobs − dcalc) ≤ ϵ) or when the difference between two consecutive iterations is small enough
(≤ ϵ/100) (Zeyen & Achauer 1997).

At the end of iterations, the resolution matrix is also computed (Zeyen & Achauer 1997; Tiberi et al. 2003). Its diagonal terms give the resolution of each separate parameter. And the off-diagonal terms quantify the inter-dependency between each parameter. This information is more complete but much more difficult to analyse because the more parameters, the biggest the matrix. That is why, we will hereafter only use the diagonal terms. Classically in this type of gravity inversion, especially with synthetic inversion of homogeneous data, the resolution is quite homogeneous inside a same layer and decrease with depth. The detail of the parameters equation can be found inZeyen & Achauer (1997); Tiberi et al. (2003).

The final density distribution mainly depends on 1) the resolving power of our method, i.e. the abil-
ity to detect small size anomalies, and 2) the assigned weight for the different constraints introduced to
stabilize the process (data and parameter variances eqs. 9 and 10, smoothing eq. 11). Moreover in the
case of GG inversions, it is necessary to take into account the complexity introduced by the simultane-
ous inversion of several, or the whole, tensor components. In theory, compared to the gravity inversion
(only Gz component), the use of the measured complete tensor in the inversion process should bring
more information on the geometry of the structures while reducing the non-uniqueness in inversion.
But in practice the behaviour and impact of the six components in the inversion have to be evaluated.
To address these issues, we test several inversions with different GG components configuration and
compare our gravity gradient inversion results with the ones from gravity (G) inversion. We first test
our approach on simple and classic synthetic tests (3-D prism geometry) to completely understand and
apprehend the subtlety of the GG inversion. This simple geometry is convenient with our parameter-
ization and creates a high frequency signature in the gradiometric signal (Okabe 1979), challenging
for the inversion. Finally we invert a more realistic synthetic case.

4 RESOLVING POWER

We first probe the accuracy of the inversion when applied to satellite data (225 km in altitude) to
define the minimal width, density contrast and maximum depth of a given body to be detected. These
informations will also guide us to define the grid inversion characteristics (cell size). Second, we
quantify the ability of the inversion to discriminate between two similar bodies by determining the
minimal distance needed to distinguish them separately.

The detection of the signal associated on a cube is mainly function of the precision of the GG
data. It only depends on the type of dataset (e.g. GOCE data only, Pail et al, 2011; or a combination
with other satellite data, Bruinsma, 2010) and on the considered tensor component (Rummel et al.
2011; Bouman et al. 2015). Due to the GOCE gradiometer configuration, the $T_{xx}$, $T_{yy}$, $T_{zz}$ and $T_{xz}$ components are more accurate (about twice better) than the $T_{xy}$ and $T_{yz}$ ones (Rummel et al. 2011). GOCE data precision being better than 10 mE, we consequently choose this value as our reference precision for the $T_{xx}$, $T_{yy}$, $T_{zz}$ and $T_{xz}$ components and half less good (20 mE) for $T_{xy}$ and $T_{yz}$. We deliberately choose a conservative precision value because it allows us to define the lower bound for the imaging capacity of GG GOCE data without being concerned about the type of data used/modeled (along the satellite track, gridded, etc.). From now on, we will consider features as detectable when their gradiometric signal exceeds those values.

## 4.1 Accuracy

We quantify the minimal characteristics (size and density contrast value) as a function of depth. From equations 3 to 8 we compute the GG response of a cube varying its size, density contrast and depth (top). For each tensor component, we follow the evolution of the extreme (maximum or minimum) value of the signal as a function of these three parameters (Fig. 2 for $T_{xx}$. See caption for more details). For this cube configuration (50 km size & 0.1 g.cm$^{-3}$), it appears that this structure is detectable on the $T_{xx}$ component (green curve) for depths above 100 km.

Figure 3 represents the evolution of the 10 mE isoline of the $T_{xx}$ component as a function of the cube dimensions and the depth of its top (see Fig. 20 for the other components). We varied the cube width between 1 and 100 km and its depth between 0 and 300 km. We also reported the curves for 4 density contrasts identified through different colors (from 0.05 up to 1 g.cm$^{-3}$). High contrasts (0.5–1 g.cm$^{-3}$) mainly correspond to crustal depth features (e.g. uncompacted sedimentary basins, variations of the Moho topography). Lower values (0.05–0.1 g.cm$^{-3}$) can be associated with lithological variations within the crust (Christensen & Mooney 1995) or thermal and/or compositional variations in the upper mantle (Afonso et al. 2013; Globig et al. 2016). We also delimited three zones; the crust (0–50 km depth), the lithospheric mantle (50–150 km depth) and the asthenosphere (150–300 km depth, Figs. 3 & 20).

For each zone and for a given density contrast, we can define the minimal size of an object to be detectable through GOCE data. For example an object with a density contrast of 0.1 g.cm$^{-3}$ located at the surface (0 km in depth) has to be at least 28 km long to be detected in the $T_{xx}$ component (Fig. 3). The same object should be wider (35 km) if located at the base of the crust (50 km in depth), and should reach ~80 km width at 300 km depth to be discerned out of noise. The required minimum cube size for each component, deduced from Fig. 20, is reported in Table 8. From these tests it appears that 1) $T_{zz}$ offers the lowest detectability threshold, 2) $T_{xx}$, $T_{xz}$ and $T_{yy}$ display a similar sensitivity, and 3) $T_{xy}$ and $T_{yz}$ display the lowest sensitivity of all components, as expected.
## 4.2 Discriminating capacity

In this section, we test the capacity of our inversion to discriminate two distinct bodies from one single large body within the GOCE signal wavelength range.

We compute the GG response for two 50 km wide cubes located at 100 km depth with a density contrast of 0.5 g.cm⁻³ and vary their separating distance (SD) along the X direction. Fig. 4.A represents the evolution of the $T_{xx}$ component along a 350 km long E-W profile (X-axis, Y=1000 km) crossing through the two bodies, as a function of the SD. For close structures (e.g. 50 km apart, dark-blue curve, Fig. 4.A), the apex of the signal is located between the two bodies and is equivalent to the signature of a single large body, as expected. As the SD increases between the two cubes (warm color curves), one can distinguish more easily the respective signature for each body, as two minimum values appear (red curves, Fig. 4.A). The emergence of those two minima defines the minimum SD value required to discriminate two bodies (noticed $SD_{req}$). We report in Fig. 4.B the evolution of the minimum value(s) of the $T_{xx}$ component (red line) and its central value along the profile (X=1000 km, black line) as a function of the SD value. For small SD value ($\le 200$ km), the signal of the two bodies is equivalent to the one of a large anomaly, and so both curves are similar. The diverging point of the two curves indicates the SD threshold value, from which we can discriminate the two bodies (noted $SD_{req}$). In our example (Fig. 4.B) $SD_{req}$ is about ~210 km. We repeated this methodology for size and depth ranges used in Fig. 3. Results are summarized in figure 5.

From those tests, we can draw several conclusions: 1) the depth of the structure (Fig. 5.B) is of higher impact on the $SD_{req}$ than its size (Fig. 5.A), and 2) for a given depth, $SD_{req}$ increases with decreasing bodies size. From Fig. 5.B we deduce a linear relationship between the $SD_{req}$ value and the depth $z$ of the body top: $SD_{min} = \frac{2}{3}z + b$, where $b$ is equal to 143 km in average and can be seen as the $SD_{req}$ for two anomalies located at the surface. As the smallest bodies give the highest $SD_{req}$ value (Fig. 5.A), we can combine this simple relationship with the minimum detectable size presented in Table. 8 to infer the minimal SD required to distinguish between two bodies. For example, 0.5 g.cm⁻³ density contrast bodies buried at 50 km depth, required a ~160 km $SD_{req}$ to discern the structures in the $T_{xx}$ signal. This value increases up to 240 km for lithospheric bodies and to 340 km for asthenospheric ones (Fig. 5.B).

An additional observation is that the density contrast value has no direct effect on the $SD_{req}$ value. Concretely it will only affect the signal amplitude. By contrast the data precision is a key parameter controlling the $SD_{req}$ value. As we can see on Fig. 4.B, a 10 mE precision for this case (density contrast of 0.5 g.cm⁻³) gives an error of around 20 km on the determination of $SD_{req}$. A precision increase directly corresponds to a $SD_{req}$ decrease (and conversely). In our study we consider a constant 10 mE error on the data regardless of the data amplitude. In this way the highest the signal amplitude,
the lowest the error. Similarly for the previous case (Fig. 4.B) if the signal amplitude increases (i.e. through the density increase of both anomalies) then this 10 mE precision will result in a lower error on the SD<sub>req</sub> determination.

# 5 SENSITIVITY OF THE INVERSION TO REGULARIZATION FACTORS

Our inversion method, as described in section 3.2, is regularized through three factors:

1. The variance of the data (hereafter VARD, $C_d$ in eq. 9) quantifies the data error distribution (Gaussian distribution). For each $T_{ij}$ component we define its variance as a percentage of the corresponding signal amplitude.

2. The variance of the model parameters (hereafter VARP, $C_p$ in eq. 10) reflects the model parameter variability. It affects the extrema values of the recovered density contrast.

3. The smoothing factor (hereafter SMO, $C_s$ in eq. 11) controls the wavelength of the structure in the model. The lower the SMO value, the smoother the model.

We test different values and combinations of VARD, VARP and SMO on the inversion of a rectangular prism of 200x200x40 km dimensions, located at 20 km depth with a density contrast of 0.2 g.cm<sup>-3</sup>. The synthesis of these inversions is represented as L-curves (Fig. 6). Each color corresponds to a different VARP value (green: 0.02, brown: 0.1, black: 0.4). High (0.001) and low (0.8) smoothing values are represented by solid and dashed lines, respectively. It is important to notice that these curves do not represent the inversion convergence, each point being a final inversion result after 5 iterations, independent of its convergence.

Global behaviours for all the components can be highlighted from Fig. 6. First a rapid decay of RMS (from ~50 to 0.1-1 mE) for a small VARD value (0-15 %) is observed. For VARD greater than ~15 % the RMS decrease is limited, except for the particular case of SMO = 0.8 & VARP = 0.4. Low smoothing (dashed lines) displays the highest values of RMS, for a given VARD value. Very small VARP values (0.02, green lines) display similar evolution, whichever the smoothing value. By excessively reducing the parameter variance, we reduce the ability of the inversion to search within the parameter space. Decreasing the smoothing constraint does not counteract this effect. The extreme opposite case, high parameter variance and low smoothing (i.e. high freedom, black dashed line, Fig. 6), leads to the highest RMS values and the presence of instabilities or no linearity in the RMS decrease with the VARD increase. For data variance (VARD) lower than ~10 %, no clear relationship can be deduced from the curves between RMS, SMO and VARP. We relate it to the simultaneous inversion of the six components which requires a certain level of freedom and a strong model coherency, not possible with these low VARD values.
A greater dispersion, in terms of RMS value, is moreover observed for $T_{xx}$, $T_{yy}$ and $T_{zz}$. This contrasting behaviour for the diagonal terms can either reflect a real difference in the sensitivity of the tensor components, or an artifact coming from the symmetry of our synthetic structure (prism), which produces symmetrical signatures.

Now, if we go into the details of these L-curves (i.e. the inversion results at each L-curve points, Fig. 6), two domains can be distinguished in terms of density contrast model, gravity gradient responses and convergence, according to the data variance value. Indeed, it appears that for VARD values lower than ~15 %, all the inversions diverge (i.e. RMS increases through iterations, Fig. 7). Even for small RMS values (<1), inversions associated with small data variance seem unlikely to retrieve a correct density model. It is only from $VARD \simeq 15\%$ that the inversion begins to correctly converge independently from the model parameter variance and the smoothing value. In those cases, density models (Fig. 8) and associated GG responses are coherent. We relate this to the nature of the gravity gradient. In this simple prism case, the associated signal is either localized on the edges ($T_{xx}$, $T_{xz}$, $T_{yz}$ and $T_{yy}$ components), the corners ($T_{xy}$) or the center of the prism ($T_{zz}$, Fig 20.f). For the two first cases (edges and corners) we clearly see that the six components can hardly be fitted simultaneously without a certain degree of freedom for the inversion. From VARD greater than 25 %, inversions converge through similar density models.

We notice that even if the gradient inversion does not totally overcome the problem of source depth location, it improves the focusing of the signal at the right depth range, compared to the gravity inversion (Fig. 8) (Mikhailov et al. 2007). This is well illustrated in the case of 2 prisms only located in the second layer of the model (20–40 km depth). The ground GG inversion results do not display any anomalies in the first layer, while the ground G inversion results present a vertical smearing in the upper layer (Fig. 8.A). Both G and GG cases displays similar vertical smearing in the lower layers. This effect is quite difficult to observe in the satellite data case (Fig. 8.B) because of the large distance between the sources and the data. The satellite data have a long wavelength content which favours a rather smooth final density model.

# 6 BENEFIT OF THE TENSOR COMPONENTS

To fully understand the benefit of the GG tensor in the inversion, a detailed study of the different $T_{ij}$ components inversion is required. However in order to get rid of the particular effect of a prism case (structure orientation, symmetry and shape of the GG signal), we applied our method on a more realistic structure. We use the Izu Bonin Marianna (IBM) slab structure derived from the slab1.0 model (Hayes et al. 2012). Such a structure is well adapted to our purpose because it displays a
large wavelength geometry coherent with the GOCE data resolution, together with lateral asymmetric variations in depth.

We model the IBM slab using the GEEC algorithm (Saraswati et al. 2019) on Hayes et al. (2012) slab geometry. This method uses irregular polyhedrons, which enables to almost perfectly reproduce the slab geometry. We subdivide the medium into several layers of different density value (according to the IASP91 model, Kennett & Engdahl 1991). We then assign for all the IBM structure and at each layer an anomaly of +5% relative to the background density value. We compute the gradiometric and gravity signatures of the slab at the same altitude (225 km) and finally invert them using the same initial homogeneous density model. We follow the procedure described above, and run gravity and gravity gradient inversions for different sets of parameterization (VARD, VARP, SMO). The gravity data variance (VARD) is classically set to 2 mGal (e.g. Tiberi et al. 2003, 2018).

A first major observation from these tests is that the inversion dependency to the VARD value is now considerably reduced. For the prism case (Figs. 6 & 8) VARD values greater than 15 % are required (section 5), whereas for IBM case, VARD of ~5–10 %, are adequate and can even lower down to ~2 % when combined with appropriate VARP and SMO values (e.g. 0.4 and 0.02, respectively, Figs. 9 & 10). It highlights the importance of the synthetic test geometry choice, which we discuss later in the article.

Three main observations stand out when comparing the gravity gradient inversion with the gravity inversion: 1) GG inversions, with appropriate VARD value, are clearly much less sensitive to the VARP and SMO values than the gravity inversion. For example, a high freedom combination of these factors (SMO = 0.8 and VARP = 0.4) leads to an unrealistic density model in the inverted gravity signal (very noisy with lots of artifacts and variations in the shape and amplitude, Fig. 11, bottom), while the GG inversion gives stable density contrasts (Fig. 11, top). 2) For the best models of both G and GG inversions (Figs. 9 & 10), independently from parameter model variance and smoothing constraint, the models derived from GG inversions always display much less artifacts and slightly better retrieve the slab geometry and the density amplitude. Finally 3) in both cases (GG and G) although the IBM model (Hayes et al. 2012) is defined until ~700 km, retrieved density anomalies are mostly located in the first 110 km, whatever the SMO, VARP and VARD values. This problem of vertical location is an inherent problem to the gravity field inversion (Li & Oldenburg 1998), since the inversion tends to raise the structures towards the surface (Fig. 8).

## 6.1 Quantification of the role of each tensor components

The simultaneous inversion of the six tensor components presents two major advantages, compared to the gravity inversion. Firstly it gives better results, hardly observable for the case of satellite data
(Figs. 8.B, 9 & 10) but undeniable for the ground data case (Fig. 8.A). Then GG inversions, when VARD value is wisely selected (section 5), are less sensitive to the other inversion parameters (SMO and VARP, Fig. 11), resulting in more robust results.

However, each component having its own accuracy (Bouman et al. 2015) and waveform characteristics, it is important to evaluate their behavior and usefulness in the GG inversion. We consider 3 cases : one component inversion (case 1), the diagonal components inversion (case 2), and the off-diagonal components inversion (case 3). At the end of the process, we also compute the response of the non-inverted components from the obtained density model (thereafter called composite response), from which we deduce a RMS (thereafter called composite RMS). Even if the RMS is first of all a mathematical criterion, in our study, it will help us to quantify the quality of a density model. In others words how a model obtained from the inversion of a part of the tensor can explain the whole tensor.

Figures 12–14 summarize all the L-curves obtained for the 3 study cases (1, 2 and 3) for the IBM geometry, and for each component. For each inversion (a point in the L-curves), we record the RMS decrease. For a given case (lines *a* to *f* in Fig. 12, Figs. 13 & 14) and for each component (rows), we calculate a mean RMS decrease (RMS<sub>dmean</sub>) from all the RMS values obtained from the inversions. All of these values are summarized in Table 2.

### 6.1.1 One component inversion (Case 1)

For a given component (Fig. 12), the highest RMS<sub>dmean</sub> always occurs when this component is inverted (for each column in Tab. 2). For example, the maximum RMS<sub>dmean</sub> for the T<sub>xx</sub> component is reached (76 %) for its own inversion (Fig. 12.a & Tab. 2). This maximum varies between 59 % (T<sub>yy</sub>) and 86 % (T<sub>xy</sub>). However this indicator is not necessary representative of the importance of a component relative to the others. Indeed several factors could be involved in a RMS variation: 1) the RMS decrease represents the difference between the initial RMS, computed from an initial model (here homogeneous) and the final RMS. If the initial model is for some reason relatively closed to the best-fitting model, then the RMS decrease will not be so high, 2) the general structure of the slab is N-S oriented (Y-axis, Fig. 9), and the highest variations of the signal will be mainly located on the derivatives of the X-component of the gravity field. This is why the highest values of average RMS<sub>dmean</sub> are observed for T<sub>xy</sub>, T<sub>xz</sub> and T<sub>xx</sub>, while the lowest value is observed for T<sub>yy</sub>. But even if this value can not be exclusively representative of a component usefulness, it can be used as a reference value (RMS<sub>ref</sub>). In other words as only one component is inverted, this percentage represents the maximal value of RMS decrease which can be reached for this component. Thus for each inversion case, the difference between RMS<sub>ref</sub> and the composite RMS<sub>dmean</sub>, noticed RMS<sub>ddif</sub>, will bring us information about the capacity for a final best-fitting model to explain all components (inverted or not). We
then use it as a proxy for the quality of the density model, even if the full tensor is not used and even if the quality of the inversion does not solely rely on a RMS value.

The detailed study of RMS<sub>dmean</sub> and RMS<sub>ddif</sub> (Tab. 2) brings essential information about the supply of each tensor component in the inversion process. If we report the mean RMS<sub>ddif</sub> of the non inverted components for each inversion case (Tab. 2.B, right row), it appears that the highest differences (worst fit) are for the inversions of the T<sub>xx</sub> and T<sub>yy</sub> components (~ 35%), followed by the T<sub>yz</sub>, T<sub>xy</sub> and T<sub>xz</sub> cases (~ 2.5 times lower) and then by T<sub>zz</sub> inversion (~5 times lower). It seems rather logical because T<sub>xx</sub> and T<sub>yy</sub> components will only bring information about structures oriented along the X- and Y-axis, respectively. Components T<sub>yz</sub>, T<sub>xy</sub> and T<sub>xz</sub> bring constrains for at least two different directions of the structure, thus the composite responses for the other components will be better than T<sub>xx</sub> and T<sub>yy</sub> inversions. Because the T<sub>zz</sub> component does not favor any orientation in the slab geometry and because it displays the highest amplitude among the tensor components, its inversion displays the lowest RMS<sub>ddif</sub> on average. T<sub>zz</sub> also displays the minimum T<sub>xx</sub> and T<sub>yy</sub> RMS<sub>ddif</sub> values (2.23 and 4.21%, respectively).

From a density distribution recovery point of view, our synthetic cases show several issues. The inversion of a specific tensor component will influence the position and orientation of artifacts that can be produced when too low VARD values are used (section 5). Logically they are aligned along the direction of the considered component used in the case of IBM (Fig. 15). The T<sub>xx</sub> and T<sub>yy</sub> component inversion create long wavelength artifacts located near the slab and oriented along the X and Y-axis, respectively (Fig. 15.a,d). With the T<sub>zz</sub> component, the artifacts are smaller and always located near the slab (Fig. 15.f). For the other components (T<sub>xy</sub>, T<sub>xz</sub>, T<sub>yz</sub>, Fig. 15.b, c, e) we observe 1) small wavelength artifacts located near the structure, and oriented along the considered components, and 2) artifacts of higher amplitude and dimensions, located on the edge of the model.

For the best models, obtained with a VARD value higher or equal to 15% clear patterns appear. In the T<sub>xx</sub> case (Fig. 16.a), the long wavelength geometry and amplitude of the IBM slab are recovered. With T<sub>yy</sub> (Fig. 16.d), a more detailed and complex geometry is modeled but the amplitude is underestimated. It seems that T<sub>zz</sub> (Fig. 16.f) is a compromise between these two characteristics with both amplitude and full geometry recovered. Similar to the RMS behaviour, the inversions of the off-diagonal components (Fig. 16.b,c,e) produce models which are a mix of the ones obtained with the corresponding components.

### 6.1.2 Diagonal components inversion (Case 2)

T<sub>xx</sub>, T<sub>yy</sub> and T<sub>zz</sub> are inverted together and the composite responses and RMS for the off-diagonal components (T<sub>xy</sub>, T<sub>xz</sub> and T<sub>yz</sub>) are computed (Tab. 2.a). First the geometry of the L-curves (Fig. 13)
for the diagonal components radically changes with an increase of the RMS with VARD. The off-diagonal L-curves display a decrease with VARD (except for $T_{yz}$). The mean $RMS_{ddif}$ is ~7%, slightly better than for the $T_{zz}$ inversion (Table 2). In detail, the lowest values are obtained for diagonal components and $T_{yz}$ (~1%). They are close to the $RMS_{ref}$ obtained by the component inversions alone. For $T_{xy}$ and $T_{xz}$ ($RMS_{ddif}$ ~ 15%), the inversion of the three diagonal components does not allow a better fit than the one obtained when $T_{zz}$ alone is considered (14 and 6%).

In terms of density distribution results, as previously observed for the prism case, low VARD values (< 15%) favor artifacts in the model (Fig. 17.a). Their dimensions and amplitude depend on VARP and SMO values. Those noisy models all produce similar good decrease of RMS, indicative of a calculated gradient signal close to the observed one. This points out the non-uniqueness of potential methods even when using several components of the gravity gradient. On the contrary best inversions (VARD > 15%, Fig. 18.a) clearly image the slab structure without artifacts near the structure. However small artifacts are present on the model boundaries.

### 6.1.3 Off-diagonal components inversion (Case 3)

As for the diagonal component, the simultaneous inversion of $T_{xy}$, $T_{xz}$ and $T_{yz}$ gives the smallest RMS for these components (respectively 6%, 2% and 0.7%, Tab. 2). It also displays a behavior similar to case 2 (an increase of the RMS with VARD). The mean $RMS_{ddif}$ (6%) is the lowest value of all cases. Similar to the previous case, the non-uniquity of the solution appears through a good decrease of RMS (Fig. 14) associated with noisy models with artifacts (Fig. 17.b). And as for case 2, best inversions (VARD > 15%, Fig. 18.b) confidently retrieve the slab structures while displaying small artifacts on the model boundaries.

The detailed study of each tensor components role in a gravity gradient inversion process highlights several major conclusions:

1. In terms of RMS value, a component is always better explained when it is inverted alone.

2. The resulting density contrast models will present slab structure and artifacts mainly oriented along the considered component.

3. The inversion of the diagonal and off-diagonal components (case 2 and 3) gives similar RMS decrease. In details, case 2 (Tab. 2) almost perfectly explains all the components except $T_{xy}$ and $T_{xz}$. In case 3, all the components are homogeneously retrieved.

4. The inversion of the 6 components (Tab. 2, last line in A and B) displays an RMS decrease slightly worse than cases 2 and 3. But this limitation is compensated by the resulting density models displaying less artifacts, particularly on the edges of the model.
5. This improvement of the inversion stability may indicate a decrease of the non-uniqueness problem. It surely highlights the importance of considering the six tensor components together in the GG inversion.

## 6.2 Bias from the synthetic case geometry

The geometry of the causative body strongly influences the way the inversion behaves. However it appears that for the gravity gradient inversion, some important bias can occur according to the synthetic case geometry. We deduced first order behaviors from the inversion of a prism structure (section 5), as it is classically proposed (e.g. Obrebski et al. 2011; Tiberi et al. 2018). However our tests on a realistic geometry brought new and additional conclusions.

When all the tensor components are involved in the inversion, a simple structure with sharp boundaries (prism) requires a high VARD value to produce a significant RMS decrease and a coherent density model (Figs. 6 & 8). A smoother and more realistic structure (IBM) displays the opposite behavior: a high VARD value leads to RMS increase (Figs. 9–14). We propose that this pattern comes from the geometry of the synthetic body. In the case of a prism, the sharp boundaries produce high frequency content in the overall gradient signal. This complexity needs an adaptative parameterization, and increasing VARD value allows to better and simultaneously fit the high frequency of each component. When the VARD value is adequate, the associated density model is coherent and corresponds well to the initial structure (Fig. 8).

The geometry of the body also influences the choice of VARP, as shown by the VARP = 0.02 case (green curves, Figs. 6 & 12–14). This value is suitable for the prism case (stable L-curves for all SMO tested values), whereas it is less appropriate for the IBM study (incoherent and unstable L-curves compared to the others VARP values). This behaviour is independent of SMO values.

Finally we showed that the choice of the synthetic case structure has multiple impacts on the inversion results either in terms of inversion convergence or final density models. First its geometry requires specific parameters values (SMO, VARP, VARD). Indeed for a prism case with sharp boundaries, inversion diverges if the VARD value is not adequate. On the contrary, the inversion of a more complex but real case (IBM slab) will display more stable inversions, besides less sensitive to the inversion parameters. However in this case the general N-S slab orientation artificially highlights some components ($T_{xx}$, $T_{xy}$, $T_{xz}$). A specific weighting of the different components will probably be necessary in the future case of real application. This weighting can be done through the choice of adequate VARD values. This weighting is strongly case dependent and needs to be thoroughly tested for each study area. It has to be adjusted in function of the data content signal, the geometry of the target and the inversion objectives. For example in the IBM case, considering the N-S orientation of the structure, it
may be interesting to give more weight to the $T_{xx}$, $T_{xy}$, $T_{xz}$ components if we want to focus on these
most "useful" components. Conversely if we want to equally explain the whole tensor, a lower weight
for these components may be required. As for our other inversion parameters (e.g. smoothing) or with
other inversion methods, the adjustment of the different parameters always requires some tests.

# 7 DISCUSSION

## 7.1 What does GOCE gravity gradient inversion bring to the Earth interior study ?

### 7.1.1 Comparison at the satellite height

As previously, we reproduced the geometry of the IBM slab from the Slab1.0 model (Hayes et al. 2012) with the GEEC algorithm (Saraswati et al. 2019) and computed the gravity and gravity gradient responses at the satellite height. We computed the gravity signal at the same height (225km) to ease the direct comparison between G and GG inversion results.

We conduct this comparison along two directions. First, we estimate the behavior of the inversions and then we compare them in terms of density contrast model recovery. For the inversion behavior, we saw that depending on the structure, the inversions requires different parameterizations, especially for VARD (section 6.2). Indeed, the inversion of the GG data at the satellite height needs a high VARD value, reflecting the ability for the inversion to depart from the data. It is especially true when the synthetic case is a prism with sharp boundaries. Too low VARD values produce diverging inversions (i.e. the RMS increase through iterations, Fig. 7) with incoherent density models and GG responses. We relate it to the nature of the gravity gradient signal, which contents small wavelengths located on the borders or corners of the body. For a smoother geometry (IBM case), this pattern is overcome and low VARD values produce coherent results (Figs. 9 & 10). Gravity better handles this VARD issue due to its signature, focused on the center of the body. The convergence of the gravity inversion does not depend on VARD. The other inversion parameter (VARP & SMO) have less influence on the GG inversion than the VARD (Fig. 11). This is probably due to the nature of the signal composed by six different signatures yet complementary (Pajot 2007). This additional information reduces the number of models which can explain data and so reduce the non-uniqueness of the inversion.

In terms of density contrasts the models issued from the GG data always better reproduce the structures compared to the gravity data only inversion. This tendency is enhanced in the simple prism inversion case with data collected at the surface, where we showed that the GG density model displays a better location of the structure and a better vertical resolution (Fig. 8.A). It appears that the use of all the tensor components overcomes an inherent problem of the poor vertical depth resolution of gravity inversions (Li & Oldenburg 1998) and its propensity to relocate sources shallower than their
real location. The combined use of the six tensor components brings more constraints to the inversion and decreases the uncertainty in the vertical location of structures. Moreover, the use of gradient can get access to the absolute density values, contrary to the gravity anomaly inversion which only gives a relative density contrast. This difference between GG and G results highlights the superiority of the GG method compared to the G one in geophysical imaging, as now used widely in applied geophysics. This feature is less straight forward for a complex and large geometry structure like IBM slab, the differences between G and GG density models remaining smaller in terms of density contrast amplitude or even recovered geometry (Figs. 9 & 10). It certainly comes from the large wavelength of the IBM structure, easily reproduced in the satellite data obtained at 225 km height. The smooth and large geometry of the IBM structure is then efficiently reproduced by both gravity and gradient data when data is recorded at satellite height.

### 7.1.2 Complementarity between satellite gravity gradient and ground gravity data

An interesting combination is to invert gravity gradient at GOCE height with ground gravity data. The GG data of lower resolution displays a global and homogeneous coverage. In spite of having generally an heterogeneous and limited data covergae, G ground data have a much higher resolution. This combination also fills the gap between local and regional structures by completing the gravitational field spectrum.

We estimate the complementarity between satellite GG data and ground G data on the realistic geometry of the IBM case. The satellite GG data are computed as previously mentioned (Saraswati et al. 2019). We build the synthetic ground gravity data set using real data coverage extracted from the BGI (Bureau Gravimetrique International) database. To limit the computation time we only consider 1/10 of the whole database, preserving the ratio between well, partially and un-covered regions. About 60 000 synthetic gravity data are computing.

Figure 19 presents the best model obtained from the inversion of these data through horizontal sections for each model layer. The important number of data generates some small size artifacts for the first layers, particularly for layer 1. As for the GG satellite model (Fig. 9), the problem of vertical location is also present. A slightly improvement of the horizontal limits of the slab is observed but is counterbalanced by the multiplication of structure with opposite density contrast sign. Overwhole G ground and GG satellite data inversions display quite similar density models (Figs. 9 & 19). Gravity gradient inversion is then an efficient tool to image lithospheric and mantle structures, while preventing field accessibility problem and new data acquisition cost especially in remote/water areas.

Nevertheless we are totally aware of several limitations in our approach. First, our test is based on a global scale slab synthetic model (Slab1.0, Hayes et al. 2012). Even if this model is deduced from
many geologic and geophysics studies all around the globe, its resolution remains poor for overturned
slabs, which is the case for IBM (Hayes et al. 2012). Moreover Slab1.0 only considers the global
and general slab structures and does not contain information about smaller and more local features
that exists in the IBM area. For example the eastern part of the Marianna trench displays a variation
of the oceanic crustal thickness ranging from 5.3 to 7 km, associated with the presence of numerous
seamounts as high as 2 or 3 km (Oakley et al. 2008). We have no clue whether these second order
structures could have been well resolved by the ground gravity data or by satellite data and their inver-
sions. Another simplification in our test is that we apply a constant density contrast value for the whole
slab. By doing so, we decrease the complexity of the structure and thus restrain the ground gravity res-
olution power in our inversion. For example in IBM the crustal thickness decreases from ~20 km in
the volcanic arc down to ~6 km at the ridges and ~10 km for the basin and the fore-arc. As previ-
ously shown and considering a 10 mE precision (section 4.1), the GOCE data are sensitive to 30 km
size structures with 0.1 g.cm⁻³ density contrast value at the surface. Consequently GG data would
not have the required resolution to image these crustal structures, contrary to the gravity ground data
which have already been proven efficient in similar contexts (e.g. Zhang et al. 2004; Montesinos et al.
2006; Welford & Hall 2007). Another limitation concerns our computation of synthetic gravity data.
We do not include data processing such as topographic effect, and corrections that could introduce
errors for different wavelengths in the ground gravity data.

Finally another limitation come from our choice to use at this time a cartesian frame which alters
the results. A first improvement in the future would be the implementation of the Tesseractal approach
(Uieda et al. 2016) in the parametrization to take into account the sphericity of the Earth. Our approach
being modular, it allows the inclusion of other informations whether it is in the form of a priori con-
straints (moho depth, density borne values, etc.) or directly included in the inversion process. This
could either concerns G or GG data at different height (e.g. ground and airborne) to complement the
gravitational field spectrum or through others methods (e.g. tomography, MT). It would allow to take
benefit of all advantages of each technique while decreasing their limits to get a more complete picture
of Earth interior. FirsIty, it will mix the very good but generally spacially limited resolution of these
methods with the global and homogeneous distribution of satellite data. Second, if for example we
coupled GG data with seismic data, the gravity gradient will constraint the absolute densities much
better than the seismic data only when the latter may better reflect the geometries.

Despite those limitations, it clearly appears that GOCE gravity gradient data are well adapted to
the study of large wavelength structures, such as Moho variation at continental scale, plates structures
and associated geodynamics). In those cases and with these objectives (geodynamic studies, structures
& geometry imaging) ground gravity measurements do not appear more indispensable than already
measured GOCE data, in comparison with their limitations (cost, distribution, etc.).

8 CONCLUSIONS

We have developed a new method to invert the gravity gradient data, based on the combination of the
forward prism computation of Uieda et al. (2016) and a Bayesian resolution approach. Because of the
high potential of the satellite gravity gradient data (GOCE mission) in the study of the Earth's interior
structure and dynamics, we chose to test our method on GOCE gravity gradient data configuration. At
this point, two main questions have guided us; which kind of target will we be able to image? And is
GOCE data inversion an interesting method for the Earth's interior study?

Considering a data precision of 10 mE and our prism parameterization, it appears that for rea-
sonable density contrast values (~0.1 g.cm⁻³) a crustal structure of around 30 km size is detectable.
Reasonably, its size needs to increase with depth for being still discernible (35 km at 50 km depth and
around 80 km at 300 km depth). The resolution can be improved through a better pre-processing of the
GG-GOCE data combined with the introduction of independent yet complementary dataset (seismic,
terrestrial / airborne gravity data).

As for the usefulness of the GOCE data inversion, we firstly quantified the benefit of each tensor
component in the inversion and estimated which combination of components is the most suitable and
then compared the inversion results with the ones we can get with the gravity inversion.

Estimating the efficiency of a given inversion is not that trivial because it depends on the factor
we consider as the most representative of the inversion quality. If we consider that the final RMS
is the most important factor, then we showed that a component is always best retrieved when solely
inverted. However if we consider that the quality of an inversion is more linked to its capacity to be the
less sensitive to the inversion parameters, then these single component configurations should not be
selected. In those cases the density models are very sensitive to the inversion parameters (SMO, VARP,
VARD) and the imaged structures are mainly oriented along the inverted components (except for T<sub>zz</sub>).
On the contrary when several components are combined together, the RMS of the inverted components
are not optimal but the ones of the composites (not inverted components) can be improved (depending
on the combination and the single component inversion). Moreover these inversions are less sensitive
to the inversion parameters, i.e. give more stable results independently from the parameters values and
consequently decrease the non-uniqueness. It highlights the interdependence of the tensor components.
Even if five of the six are mathematically independent (the trace tensor being null, T<sub>xx</sub>, T<sub>yy</sub> and T<sub>zz</sub>
are dependent), when used in an inversion to characterize a common objet each component is linked
to the others (e.g. T<sub>xy</sub> inverted value will affected T<sub>xx</sub>, T<sub>yy</sub>, T<sub>zz</sub> and T<sub>yz</sub> ones).
This conclusion has however to be weighted for the VARD dependence particularly when the method is applied on synthetic prism configuration. A too low value (ie. a small possibility to the inversion to be far from the data) leads to diverging inversions (RMS increases through iterations) with incoherent density models and associated responses. We make the hypothesis that this specifically comes from the nature of the GG signal, a same object being characterized by six different signatures. For the problematic prism case, five over the six components (except Tzz) display a high frequency signal located on the sides or vertices of the structure. A rather high VARD value counterbalances this effect and makes the inversions converge.

Finally we highlighted the interest of the GG inversion with and compared to the gravity inversion. Indeed at the same altitude, the inversion of GG data, with our method, better recovers the structures and displays less artifacts together with a better vertical location, one of the major problems with the gravity field inversion. This improvement is more significant in the case of ground data (Fig. 8.A) than satellite configuration data (Fig. 8.B). In the case of large wavelength geodynamics studies (e.g. subduction) where small and local structures are not targeted, the use of GG satellite data seems totally adapted and does not seem to bring less information than gravity ground data alone.

## ACKNOWLEDGEMENTS

This study was supported by the "Centre National d'Etudes Spatiales" (CNES) with a post-doctoral grant and a TOSCA project financing (PIGGS project). We are grateful for the enriching discussions with Rodolphe Cattin and Stephane Mazzotti. We also thank the BGI (Bureau Gravimetrique International) for the gravity data. We want to deeply thank the editor Professor Bert Vermeersen, Martin Fuchs and one anonymous reviewer for their constructive remarks, which have greatly improved the manuscript. This papier is IPGP Contribution 4147.
REFERENCES

Afonso, J., Fullea, J., Griffin, W., Yang, Y., Jones, A., D. Connolly, J., & O'Reilly, S., 2013. 3-d multiobserv-
able probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle.
i: A priori petrological information and geophysical observables, *Journal of Geophysical Research: Solid
Earth*, **118**(5), 2586–2617.

Afonso, J. C., Salajegheh, F., Szwillus, W., Ebbing, J., & Gaina, C., 2019. A global reference model of the lithosphere and upper mantle from joint inversion and analysis of multiple data sets, *Geophysical Journal International*, **217**(3), 1602–1628.

Basuyau, C., Diament, M., Tiberi, C., Hetenyi, G., Vergne, J., & Peyrefitte, A., 2013. Joint inversion of teleseismic and goce gravity data: application to the himalayas, *Geophysical Journal International*, **193**(1), 149–160.

Beiki, M., 2010. Analytic signals of gravity gradient tensor and their application to estimate source location, *Geophysics*, **75**(6), 159–174.

Beiki, M. & Pedersen, L. B., 2011. Window constrained inversion of gravity gradient tensor data using dike and contact models, *Geophysics*, **76**(6), 159–172.

Bell, R. E., Anderson, R., & Pratson, L., 1997. Gravity gradiometry resurfaces, *The Leading Edge*, **16**(1), 55–59.

Bouman, J., Rispens, S., Gruber, T., Koop, R., Schrama, E., Visser, P., Tscherning, C. C., & Veicherts, M., 2009. Preprocessing of gravity gradients at the goce high-level processing facility, *Journal of Geodesy*, **83**(7), 659–678.

Bouman, J., Fiorot, S., Fuchs, M., Gruber, T., Schrama, E., Tscherning, C., Veicherts, M., & Visser, P., 2011. Goce gravitational gradients along the orbit, *Journal of Geodesy*, **85**(11), 791.

Bouman, J., Ebbing, J., Meekes, S., Fattah, R. A., Fuchs, M., Gradmann, S., Haagmans, R., Lieb, V., Schmidt, M., Dettmering, D., et al., 2015. Goce gravity gradient data for lithospheric modeling, *International Journal of Applied Earth Observation and Geoinformation*, **35**, 16–30.

Bouman, J., Ebbing, J., Fuchs, M., Sebera, J., Lieb, V., Szwillus, W., Haagmans, R., & Novak, P., 2016. Satellite gravity gradient grids for geophysics, *Scientific reports*, **6**, 21050.

Braitenberg, C., 2015. Exploration of tectonic structures with goce in ari ca and across-continents, *International Journal of Applied Earth Observation and Geoinformation*, **35**, 88–95.

Braitenberg, C., Mariani, P., Ebbing, J., & Sprlak, M., 2011. The enigmatic chad lineament revisited with global gravity and gravity-gradient fields, *Geological Society, London, Special Publications*, **357**(1), 329–341.

Bruinsma, S., 2010. Goce gravity field recovery by means of the direct numerical method, in *ESA Living Planet Symposium 2010, Bergen, Norway, June 27-July 2*.

Cadio, C., Panet, I., Davaille, A., Diament, M., Métivier, L., & de Viron, O., 2011. Pacific geoid anomalies revisited in light of thermochemical oscillating domes in the lower mantle, *Earth and Planetary Science Letters*, **306**(1-2), 123–135.
Cambiotti, G. & Sabadini, R., 2012. A source model for the great 2011 tohoku earthquake (mw= 9.1) from inversion of grace gravity data, *Earth and Planetary Science Letters*, **335**, 72–79.

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

Dogru, F., Pamukcu, O., Gonenc, T., & Yildiz, H., 2018. Lithospheric structure of western anatolia and the aegean sea using goce-based gravity field models., *Bollettino di Geofisica Teorica ed Applicata*, **59**(2).

Douch, K., Christophe, B., Foulon, B., Panet, I., Pajot-Métivier, G., & Diament, M., 2014. Ultra-sensitive electrostatic planar acceleration gradiometer for airborne geophysical surveys, *Measurement Science and Technology*, **25**(10), 105902.

Dransfield, M. & Milkereit, B., 2007. Airborne gravity gradiometry in the search for mineral deposits, in *Proceedings of exploration*, vol. 7, pp. 341–354.

Drinkwater, M. R., Haagmans, R., Muzi, D., Popescu, A., Floberghagen, R., Kern, M., & Fehringer, M., 2006. The goce gravity mission: Esa's first core earth explorer, in *Proceedings of the 3rd international GOCE user workshop*, pp. 6–8, Citeseer.

Forsberg, R. & Olesen, A. V., 2010. Airborne gravity field determination, in *Sciences of geodesy-I*, pp. 83–104, Springer.

Fuchs, M. J. & Bouman, J., 2011. Rotation of goce gravity gradients to local frames, *Geophysical Journal International*, **187**(2), 743–753.

Geng, M., Huang, D., Yang, Q., & Liu, Y., 2014. 3d inversion of airborne gravity-gradiometry data using cokriging, *Geophysics*, **79**(4), G37–G47.

Globig, J., Fernández, M., Torne, M., Vergés, J., Robert, A., & Faccenna, C., 2016. New insights into the crust and lithospheric mantle structure of Africa from elevation, geoid, and thermal analysis, *Journal of Geophysical Research: Solid Earth*, **121**(7), 5389–5424.

Hassan, H. H., 2005. Empirical mode decomposition (emd) of potential field data: airborne gravity data as an example, in *SEG Technical Program Expanded Abstracts 2005*, pp. 704–706, Society of Exploration Geophysicists.

Hayes, G. P., Wald, D. J., & Johnson, R. L., 2012. Slab1.0: A three-dimensional model of global subduction zone geometries, *Journal of Geophysical Research: Solid Earth*, **117**(B1).

Hirt, C., 2014. Goce’s view below the ice of antarctica: Satellite gravimetry confirms improvements in bedmap2 bedrock knowledge, *Geophysical Research Letters*, **41**(14), 5021–5028.

Jacoby, W. & Smilde, P. L., 2009. *Gravity interpretation: fundamentals and application of gravity inversion and geological interpretation*, Springer Science & Business Media.

Janjić, T., Schröter, J., Savcenko, R., Bosch, W., Albertella, A., Rummel, R., & Klatt, O., 2012. Impact of combining grace and goce gravity data on ocean circulation estimates, *Ocean Science*, **8**(1), 65–79.

Kennett, B. & Engdahl, E., 1991. Traveltimes for global earthquake location and phase identification, *Geophysical Journal International*, **105**(2), 429–465.

Knudsen, P., Bingham, R., Andersen, O., & Rio, M.-H., 2011. A global mean dynamic topography and ocean
circulation estimation using a preliminary goce gravity model, *Journal of Geodesy*, **85**(11), 861-879.

Li, H., Fang, J., & Braitenberg, C., 2017. Lithosphere density structure beneath the eastern margin of the tibetan plateau and its surrounding areas derived from goce gradients data, *Geodesy and Geodynamics*, **8**(3), 147-154.

Li, Y., 2001. Processing gravity gradiometer data using an equivalent source technique, in *SEG Technical Program Expanded Abstracts 2001*, pp. 1466-1469, Society of Exploration Geophysicists.

Li, Y. & Oldenburg, D. W., 1998. 3-d inversion of gravity data, *Geophysics*, **63**(1), 109-119.

Martinez, C., Li, Y., Krahenbuhl, R., & Braga, M. A., 2012. 3d inversion of airborne gravity gradiometry data in mineral exploration: A case study in the quadrilátero ferrífero, brazil, *Geophysics*, **78**(1), B1-B11.

Menke, W., 1984. *Geophysical data analysis: Discrete inverse theory*, Academic press.

Migliaccio, F., Reguzzoni, M., Sanso, F., Tscherning, C., Veicherts, M., et al., 2010. Goce data analysis: the space-wise approach and the first space-wise gravity field model, in *Proceedings of the ESA living planet symposium*, vol. 28, Citeseer.

Mikhailov, V., Pajot, G., Diament, M., & Price, A., 2007. Tensor deconvolution: A method to locate equivalent sources from full tensor gravity data, *Geophysics*, **72**(5), I61-I69.

Montesinos, F., Arnoso, J., Benavent, M., & Vieira, R., 2006. The crustal structure of el hierro (canary islands) from 3-d gravity inversion, *Journal of Volcanology and Geothermal Research*, **150**(1-3), 283-299.

Nagy, D., Papp, G., & Benedek, J., 2000. The gravitational potential and its derivatives for the prism, *Journal of Geodesy*, **74**(7-8), 552-560.

Oakley, A., Taylor, B., & Moore, G., 2008. Pacific plate subduction beneath the central mariana and izu-bonin fore arcs: New insights from an old margin, *Geochemistry, Geophysics, Geosystems*, **9**(6).

Obrebski, M., Allen, R. M., Pollitz, F., & Hung, S.-H., 2011. Lithosphere-asthenosphere interaction beneath the western united states from the joint inversion of body-wave traveltimes and surface-wave phase velocities, *Geophysical Journal International*, **185**(2), 1003-1021.

Okabe, M., 1979. Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies, *Geophysics*, **44**(4), 730-741.

O'Donnell, J., Adams, A., Nyblade, A., Mulibo, G., & Tugume, F., 2013. The uppermost mantle shear wave velocity structure of eastern Africa from rayleigh wave tomography: constraints on rift evolution, *Geophysical Journal International*, **194**(2), 961-978.

Pail, R., Bruinsma, S., Migliaccio, F., Förste, C., Goiginger, H., Schuh, W.-D., Höck, E., Reguzzoni, M., Brockmann, J. M., Abrikosov, O., et al., 2011. First goce gravity field models derived by three different approaches, *Journal of Geodesy*, **85**(11), 819.

Pajot, G., 2007. *Caractérisation, analyse et interprétation des données de gradiométrie en gravimétrie*, Ph.D. thesis, Institut de physique du globe de paris-IPGP.

Pajot, G., De Viron, O., Diament, M., Lequentrec-Lalancette, M.-F., & Mikhailov, V., 2008. Noise reduction through joint processing of gravity and gravity gradient data, *Geophysics*, **73**(3), I23-I34.

Panet, I., Pajot-Métivier, G., Greff-Lefftz, M., Métivier, L., Diament, M., & Mandea, M., 2014. Mapping the
mass distribution of earth's mantle using satellite-derived gravity gradients, *Nature Geoscience*, **7**(2), 131.

Pedersen, L. B. & Rasmussen, T. M., 1990. The gradient tensor of potential field anomalies: Some implications on data collection and data processing of maps, *Geophysics*, **55**(12), 1558–1566.

Reguzzoni, M. & Sampietro, D., 2012. Moho estimation using goce data: a numerical simulation, in *Geodesy for Planet Earth*, pp. 205–214, Springer.

Rummel, R., Yi, W., & Stummer, C., 2011. Goce gravitational gradiometry, *Journal of Geodesy*, **85**(11), 777.

Sampietro, D., Reguzzoni, M., & Braitenberg, C., 2014. The goce estimated moho beneath the tibetan plateau and himalaya, in *Earth on the Edge: Science for a Sustainable Planet*, pp. 391–397, Springer.

Saraswati, A. T., Cattin, R., Mazzotti, S., & Cadio, C., 2019. New analytical solution and associated software for computing full-tensor gravitational field due to irregularly shaped bodies, *Journal of Geodesy*, pp. 1–17.

Tian, Y., Ke, X., & Wang, Y., 2019. A folding calculation method based on the preconditioned conjugate gradient inversion algorithm of gravity gradient tensor, *Pure and Applied Geophysics*, **176**(1), 215–234.

Tiberi, C., Diament, M., Déverchère, J., Petit-Mariani, C., Mikhailov, V., Tikhotsky, S., & Achauer, U., 2003. Deep structure of the baikal rift zone revealed by joint inversion of gravity and seismology, *Journal of Geophysical Research: Solid Earth*, **108**(B3).

Tiberi, C., Gautier, S., Ebinger, C., Roecker, S., Plasman, M., Albaric, J., Déverchère, J., Peyrat, S., Perrot, J., Wambura, R. F., et al., 2018. Lithospheric modification by extension and magmatism at the craton-orogenic boundary: North tanzania divergence, east Africa, *Geophysical Journal International*, **216**(3), 1693–1710.

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

van der Meijde, M., Pail, R., Bingham, R., & Floberghagen, R., 2015a. Goce data, models, and applications: A review, *International Journal of Applied Earth Observation and Geoinformation*, **35**, 4–15.

van der Meijde, M., Pail, R., Bingham, R., & Floberghagen, R., 2015b. Goce data, models, and applications: A review, *International Journal of Applied Earth Observation and Geoinformation*, **35**, 4–15.

Welford, J. K. & Hall, J., 2007. Crustal structure of the newfoundland rifted continental margin from constrained 3-d gravity inversion, *Geophysical Journal International*, **171**(2), 890–908.

Zeyen, H. & Achauer, U., 1997. Joint inversion of teleseismic delay times and gravity anomaly data for regional structures, in *Upper Mantle Heterogeneities from active and passive seismology*, pp. 155–168, Springer.

Zhang, J., Wang, C.-Y., Shi, Y., Cai, Y., Chi, W.-C., Dreger, D., Cheng, W.-B., & Yuan, Y.-H., 2004. Three-dimensional crustal structure in central Taiwan from gravity inversion with a parallel genetic algorithm, *Geophysics*, **69**(4), 917–924.
<table><caption>**Table 1.** Summary of the minimal cube size detectable (in km) in each GOCE tensor component, for three distinct depth ranges (0–50, 50–150 and 150–300 km) and four density contrast values (0.05 g.cm⁻³: grey; 0.1 g.cm⁻³: blue; 0.5 g.cm⁻³: red and 1.0 g.cm⁻³: black) and for each tensor components.</caption><thead><tr><th>Tensor component</th><th>Crust (0-50 km)</th><th>Lithospheric mantle (50-150 km)</th><th>Asthenosphere (150-300 km)</th></tr></thead><tbody><tr><td rowspan="4">T<sub>XX</sub></td><td>35 - 45</td><td>45 - 65</td><td>65 - &gt; 100</td></tr><tr><td>28 - 35</td><td>35 - 50</td><td>50 - 80</td></tr><tr><td>15 - 20</td><td>20 - 28</td><td>28 - 45</td></tr><tr><td>13 - 15</td><td>15 - 22</td><td>22 - 35</td></tr><tr><td rowspan="4">T<sub>XY</sub></td><td>79 - 90</td><td>90 - &gt; 100</td><td>&gt; 100</td></tr><tr><td>55 - 69</td><td>69 - 93</td><td>93 - &gt; 100</td></tr><tr><td>30 - 38</td><td>38 - 52</td><td>52 - 75</td></tr><tr><td>24 - 30</td><td>30 - 40</td><td>40 - 57</td></tr><tr><td rowspan="4">T<sub>XZ</sub></td><td>37 - 45</td><td>45 - 62</td><td>62 - 86</td></tr><tr><td>29 - 35</td><td>35 - 49</td><td>49 - 68</td></tr><tr><td>16 - 20</td><td>20 - 26</td><td>26 - 39</td></tr><tr><td>13 - 15</td><td>15 - 21</td><td>21 - 30</td></tr><tr><td rowspan="4">T<sub>YY</sub></td><td>35 - 44</td><td>44 - 61</td><td>61 - 93</td></tr><tr><td>28 - 35</td><td>35 - 48</td><td>48 - 72</td></tr><tr><td>15 - 20</td><td>20 - 27</td><td>27 - 41</td></tr><tr><td>13 - 15</td><td>15 - 21</td><td>21 - 33</td></tr><tr><td rowspan="4">T<sub>YZ</sub></td><td>48 - 59</td><td>59 - 80</td><td>80 - &gt; 100</td></tr><tr><td>37 - 45</td><td>45 - 61</td><td>61 - 87</td></tr><tr><td>21 - 26</td><td>26 - 35</td><td>35 - 50</td></tr><tr><td>16 - 20</td><td>20 - 28</td><td>28 - 39</td></tr><tr><td rowspan="4">T<sub>ZZ</sub></td><td>28 - 35</td><td>35 - 50</td><td>50 - 75</td></tr><tr><td>22 - 28</td><td>28 - 38</td><td>38 - 59</td></tr><tr><td>13 - 15</td><td>15 - 21</td><td>21 - 34</td></tr><tr><td>10 - 12</td><td>12 - 16</td><td>16 - 26</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>**Table 2.** Summary of the all A) RMS<sub>dmean</sub> and B) RMS<sub>ddif</sub>, computed for cases 1, 2 and 3 and for each component. RMS<sub>dmean</sub> is computed from all the RMS values obtained from the whole inversions realized. RMS<sub>ddif</sub> corresponds to the difference between RMS<sub>ref</sub> and the composite RMS<sub>dmean</sub>. See text (section 6.1) for explanations.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Explicative scheme for the notations used in the calculation of the gravity gradient (eqs. 2–8).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** **A**: $T_{xx}$ component response for a cube of 50 km size (white box), at 0 km depth and with a density contrast of 0.1 g.cm⁻³. The white square represents the cube location and the black line is the E-W profile displayed in **B**. **B**: E-W section in **A**. Each color line represents $T_{xx}$ for a given depth of the cube top (in km). The grey area corresponds to the undetectable signal area ie where the signal amplitude is lower than 10 mE (see text for explanations).</figcaption></figure>
<figure><img src="image_5.png" /></figure>

**Figure 3.** Position and evolution of the contour line 10 mE as a function of the cube dimensions (Y axis) and depth of its top (X axis) for the $T_{xx}$ component (see section 4.1 for more explanations). Each color curve corresponds to a different density contrast value (0.05 g.cm$^{-3}$: grey; 0.1 g.cm$^{-3}$: blue; 0.5 g.cm$^{-3}$: red and 1.0 g.cm$^{-3}$: black). The surface below a given curve delimites undetectable cases (signal < 10 mE)
<figure><img src="image_6.png" /><figcaption>**Figure 4.** **A:** E-W section of $T_{xx}$ component for the case of two identical cubes (50 km wide, at 100 km depth with a density contrast of $0.5 \text{ g.cm}^{-3}$), aligned along the X direction. Each color curve corresponds to a given distance between these two cubes (SD). The two framed insets correspond to the $T_{xx}$ component for SD = 50 km (blue) and 400 km (brown). **B:** Evolution of the minimum of the $T_{xx}$ component (black) and its center value along the profile at X = 1000 km (red) as a function of the distance between the two anomalies. When the center value (red curve) is greater than the minimum of the signal (black curve), we distinctly observe the signature of each anomaly. The red area indicates the GOCE precision (see text for explanations).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** SD<sub>req</sub> required to separately observe two distinct anomalies for all sizes and depth ranges used in Fig. 3. SD as a function of A) the cube size with color corresponding to the body depth and B) the body depth with color corresponding to the cube size.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** L-curves obtained for each tensor component from the inversion of the signature of a simple 200*200*40 km prism, located at 20 km depth with a density contrast of 0.2 g.cm<sup>-3</sup>. Each curve point represents the final inversion RMS after 5 iterations.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

**Figure 7.** Importance of the VARD value through the plots of the inversion divergence (a), final density model (b) and its gravity gradient response (c) for the case of a cube (section 5) and for VARD = 4 %, corresponding to a diverging inversion case. The result of an inversion made with an adequate VARD value (>15 %) is shown in Fig 8.B.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Inversion results of two prisms with opposite density contrasts for ground (A) and satellite (B) gravity and gravity gradient synthetic data.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Horizontal section of the density contrast model for GG (top) and gravity (bottom) inversion for the two inversion parameter sets (VARP = 0.02, VARD = 20% for GGI and 2 mGal for GI, SMO = 0.8) giving the best models.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Vertical EW sections of the two models in Fig 9 along the axis y = 3500 for the GG (top) and the gravity (bottom) inversion.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Exemple of gravity gradient (top) and graviy (bottom) density models obtained from the inversion of the IBM slab signature with a SMO = 0.8 and VAR = 0.4.</figcaption></figure>
Case 1

<figure><img src="image_14.png" /><figcaption>**Figure 12.** The whole L-curves computed in the IBM inversion, for each inversion of the case 1 (horizontal) and for each tensor component (vertical). The more visible plots shows the inverted component while the more transparent ones corresponds to the composite L-curves (see text for explanations). The color of the line corresponds to different VARP values (green : 0.02, red : 0.1, black : 0.4). Dashed and solid lines correspond to SMO value of 0.8 and 0.001, respectively.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** The whole L-curves computed in the IBM inversion, for the case 2 (horizontal) and for each tensor component (vertical). The more visible plots show the inverted components while the more transparent ones correspond to the composite L-curves (see text for explanations). The color of the line corresponds to different VARP values (green : 0.02, red : 0.1, black : 0.4). Dashed and solid lines correspond to SMO value of 0.8 and 0.001, respectively.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** The whole L-curves computed in the IBM inversion, for the case 3 (horizontal) and for each tensor component (vertical). The more visible plots show the inverted components while the more transparent ones correspond to the composite L-curves (see text for explanations). The color of the line corresponds to different VARP values (green : 0.02, red : 0.1, black : 0.4). Dashed and solid lines correspond to SMO value of 0.8 and 0.001, respectively.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Example of noisy density contrast models (low VARD values) obtained for all the GG inversions of the case 1 for the IBM inversion</figcaption></figure>
Case 1

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Example of best density contrast models (VARD value ≥ 15 %) obtained for all the GG inversions of the case 1 for the IBM inversion</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Example of noisy density contrast models (low VARD values) obtained for all the GG inversions of the cases 2 (a) & 3 (b) for the IBM inversion</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Example of best density contrast models (VARD value ≥ 15 %) obtained for all the GG inversions of the cases 2 (a) & 3 (b) for the IBM inversion</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Horizontal section of the density contrast model obtained from the inversion of the BGI database gravity data. In black the position of the slab. In white (on the right) position of the gravity data.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 20.** SUPPLEMENTARY Same as Fig. 3 but for the others tensor components. In the bottom right corner, example of the gravity gradient signature for a simple cube of 50 km dimensions, located at 50 km depth and with a density contrast of 0.5.</figcaption></figure>