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

The emergence of high resolution satellite measurements of the gravitational ﬁeld (GOCE mission) offers promising perspectives for the study of the Earth’s interior. These new data call for the development of innovant analysis and interpretation methods. Here we combine a forward prism computation with a Bayesian resolution approach to invert for these gravity gradient data conﬁguration. We apply and test our new method on satellite data conﬁguration, i.e. 225 km height with a global and homogeneous geographic distribution. We ﬁrst 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 − 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 ﬁve 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 (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 inversion parameterization (smoothing, variances, etc.). On the contrary, the simultaneous inversion of several components displays a signiﬁcant improvement for the global tensor 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 succesfully 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 nonexistance 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.


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 innovant 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 −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 (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 inversion parameterization (smoothing, variances, etc.). On the contrary, the simultaneous inversion of several components displays a significant improvement for the global tensor 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 succesfully 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 nonexistance 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

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
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 heigh 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 − → 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 important 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 resolution 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.

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. Classicaly, people prefer to use the Local North-Oriented Frame (LNOF) (Bouman et al. 2009), 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 Inversion of the GOCE gravity gradients 5 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 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 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.

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 observated and predicted GG data.

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

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 Inversion of the GOCE gravity gradients 7 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 : 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) : 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).

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: and for the parameter : The symbol T denotes the transpose operation, p 0 is the density starting model and p T = (∆ρ) 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: Equation 11 relates adjacent density contrasts ∆p with their corresponding distance ∆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: 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 = i): 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 ( ) is Inversion of the GOCE gravity gradients 9 reached ((d obs −d calc ) ≤ ) 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. Classicaly 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 ability 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 simultaneous inversion of several, or the whole, tensor components. In theory, compared to the gravity inversion (only G z 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 parameterization and creates a high frequency signature in the gradiometric signal (Okabe 1979), challenging for the inversion. Finally we invert a more realistic synthetic case.

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.

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

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 −3 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, darkblue 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 (≤ 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 = 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 −3 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 controling 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 −3 ) 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. throught the density increase of both anomalies) then this 10 mE precision will result in a lower error on the SD req determination.

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.

Inversion of the GOCE gravity gradients 13
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 theses 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 V ARD 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.

BENEFIT OF THE TENSOR COMPONENTS
To fully understand the benefit of the GG tensor in the inversion, a detailled 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. 2003Tiberi et al. , 2018 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).

Quantification of the role of each tensor components
The simultaneous inversion of the six tensor components presents two majors advantages, compared to the gravity inversion. Firstly it gives better results, hardly observable for the case of satellite data Inversion of the GOCE gravity gradients 15 (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 offdiagonal 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. All of these values are summarized in Table 2.

One component inversion (Case 1)
For a given component (Fig. 12), the highest RMS dmean always occurs when this component is inverted (for each column in Tab. 2). For example, the maximum RMS dmean for the T xx component is reached (76 %) for its own inversion ( Fig. 12.a & Tab. 2). This maximum varies between 59 % (T yy ) and 86 % (T xy ). 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 deriva- The detailled study of RMS dmean and RMS ddif (Tab. 2) brings essential information about the supply of each tensor component in the inversion process. If we report the mean RMS ddif 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 xx and T yy components (∼ 35%), followed by the T yz , T xy and T xz cases (∼ 2.5 times lower) and then by T zz inversion (∼5 times lower). It seems rather logical because T xx and T yy components will only bring information about structures oriented along the X-and Y-axis, respectively. Components T yz , T xy and T xz bring constrains for at least two different directions of the structure, thus the composite responses for the other components will be better than T xx and T yy inversions. Because the T zz 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 ddif on average. T zz also displays the minimum T xx and T yy RMS ddif 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 xx and T yy 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 zz component, the artifacts are smaller and always located near the slab (Fig. 15.f). For the other components (T xy , T xz , T yz , 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 xx case (Fig. 16.a), the long wavelength geometry and amplitude of the IBM slab are recovered. With T yy (Fig. 16.d), a more detailled and complex geometry is modeled but the amplitude is underestimated. It seems that T zz (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 offdiagonal components ( Fig. 16.b,c,e) produce models which are a mix of the ones obtained with the corresponding components.

Diagonal components inversion (Case 2)
T xx , T yy and T zz are inverted together and the composite responses and RMS for the off-diagonal components (T xy , T xz and T yz ) are computed (Tab. 2.a). First the geometry of the L-curves (Fig. 13) Inversion of the GOCE gravity gradients 17 for the diagonal components radically changes with an increase of the RMS with VARD. The offdiagonal 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.

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-unicity 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 detailled 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.

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.

This improvement of the inversion stability may indicate a decrease of the non-uniqueness prob-
lem. It surely highlights the importance of considering the six tensor components together in the GG inversion.

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 . 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,. 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 throught the choice of adequate VARD values. This weighting is strongly case dependent and needs to be thouroughly 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 Inversion of the GOCE gravity gradients 19 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 equaly 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 adjustement of the different parameters always requires some tests.

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

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.

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 inversions. 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 resolution 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 previously 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 −3 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 Tesseroidal 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 constraints (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. Firslty, 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.).

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 reasonable density contrast values (∼0.1 g.cm −3 ) 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 zz ).
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 xx , T yy and T zz are dependent), when used in an inversion to characterize a common objet each component is linked to the others (e.g. T xy inverted value will affected T xx , T yy , T zz and T yz 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.    with a density contrast of 0.5 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 distincly observe the signature of each anomaly. The red area indicates the GOCE precision (see text for explanations).