“ DIPPING PRISM MODELLING OF SUBDUCTION PLATES IN VIEW OF AN IMPROVED GOCE GLOBAL MOHO: THE TONGA EXAMPLE „

The study of subduction zones, i.e. the process occurring at convergent boundaries by which


INTRODUCTION
The study of the main subduction plates is usually performed by combining several different methodolo− gies and observations: high accuracy mapping of seis− micity [Engdahl et al., 1998;Deal et al., 1999, and updates] and moment tensor catalogues (gCMT, http://www.globalcmt.org) can, for instance, give direct information about the Benioff−Wadati zone (i.e. the boundary between the subducting oceanic lithosphere beneath either a continental lithosphere or another oceanic lithosphere), while active−source seismic data allow constraining the shallow section of the subduct− ing slabs in a similar way. The incorporation of deep seismicity (where all earthquakes occur within sub− ducted slabs, regardless of the mechanism) offers the possibility to improve the modelling by introducing a possible curvature of subducting slabs occurring over their descend into the mantle [England et al., 2004;Syracuse and Abers, 2006]. Finally, the availability of high−resolution bathymetry and sediment layers can help in defining the geometry of the slab at the trenches. One of the most complete products generated by these kinds of studies is the Slab1.0 model [Hayes et al., 2012], which is a global model that describes the three−dimensional geometry of the main subduction zones obtained by compiling and exploiting many sub− duction−related seismic data.
Since there is, generally, a large mass density dis− continuity in correspondence with subducting plates, due to the presence of a net separation between the lighter crustal and the heavier mantle masses, gravity data can be used to study subduction zones. For in− stance, ground gravity data, i.e. Bouguer anomalies, were used to study the crustal thickness and the sub− ducting lithosphere in Greece showing a good agree− ment with the recent seismological results [Tsokas and Hansen, 1997], while gravity field models derived from the GEOdetic SATellite (GEOSAT), SEAfaring SATellite (SEASAT), and European Remote−Sensing (ERS−1) satellite altimetry data were combined with multi− channel seismic records and shipborne gravity data to study the South East Pacific−Antarctic region, leading to a model of converging oceanic crustal blocks with a partial subduction in the Bellingshausen Sea [Gohl et al., 1997].
The use of gravity observations in such a kind of studies has been up to now quite limited due to the spe− cific geometry of the problem. In fact, subduction plates usually extend in large areas, thus requiring long and expensive airborne gravity campaigns. Even more, the significant depth of the slabs, which can reach several hundreds of kilometres, makes the characterization of subducting plates a difficult task. With the advent of the European Space Agency (ESA) mission GOCE (Grav− ity field and steady−state Ocean Circulation Explorer) [Drinkwater et al., 2003], new important information can be added to improve the problem solution. Actually, several studies have already shown the capabilities of GOCE data to retrieve Solid Earth information, see for instance [Bagherbandi et al., 2013;Sampietro et al., 2014;Eshagh et al., 2016;Barzaghi et al., 2016;Ye et al., 2016;Abrehdary et al., 2017;Tondi et al., 2017;Sampi− etro et al., 2018−b;Sobh et al., 2018−b;Reguzzoni et al., 2019] and in particular to infer the mantle mass density distribution. For the lat− ter topic, see Bouman et al. [2015], where GOCE grav− ity gradients have been used to determine the structure of the crust and the upper mantle in the North East At− lantic region and in Saudi Arabia, or Panet et al. [2014], where global anomaly grids of the Earth's gravitational gradients at satellite altitude have been used to study the geometry of the mantle masses down to mid−man− tle depths.
While looking at GOCE grids of gravity gradients at satellite altitude, a very high−frequency spatial compo− nent characterized by a large amplitude can be observed on the diagonal components of the Marussi tensor over subduction zones, thus leading from the one hand to the inquiry about the capability of using GOCE to study and map these regions (see Figure 1), and on the other hand to the necessity of modelling these regions, even in a simplified way, when studying the global behavior of the Earth crust from gravity observations. Here, we re− call that the Marussi tensor is the gravity gradient ten− sor, which is traceless and symmetric. In the present research, subducting plates are modelled by dipping prisms. Because of its simplicity, this modelling can be used on a global basis and will give us the chance to improve the GEMMA global Moho  by better reducing the GOCE data be− fore inverting them. According to the dipping prism modelling, GOCE gravity gradients along the orbit are used to estimate the shape, size, and orientation param− eters of the studied subduction plate along with its den− sity contrast. Note that, even if in correspondence of the subducting plates several density contrasts are present, e.g. the one related to the accretionary wedge, or to the magma melting and diapirism [Lallemand et al., 1992;Behn et al., 2011], here we will consider just the con− trast between the subducting crust and the surrounding lithosphere (see Figure 2). This approximation is justified by the fact that a detailed modelling of each subduction zone for global studies is a very difficult task, which is in general not required due to the low spatial resolution of such global crustal studies.
The proposed algorithm is tested by studying the Tonga subduction zone, located where the Pacific plate subducts beneath the Australian one at the Tonga and Kermadec trenches and dips to the West direction gen− erating one of the most active tectonic complexes in the world [Van Der Hilst, 1995]. The Kermadec−Tonga sub− duction zone is characterized as a convergent plate boundary, which stretches from the North Island of New Zealand northward, and includes the Hikurangi Trough, and the Kermadec and Tonga Trenches [Ewart et al., 1977]. Along this zone, the Pacific plate to the East is subducting beneath the Indo−Australian plate at a rate of 5.5 to 7.4 centimetres per year [Garcia−Castellanos et al., 2000].
In Section 2, the data reduction to isolate the sub− duction anomalous signal and its inversion using a simulated annealing (SA) procedure are presented. In Section 3, the developed algorithm is applied on syn− thetic data in order to assess its performance in recov− ering the unknown subduction parameters and to eval− uate the contribution of every single component of the gravity gradient tensor within the inversion procedure. Section 4 addresses the real−data processing and the consistency of the obtained results with those coming from existing seismic models. A final discussion on the obtained results is held in Section 5, confirming the feasibility of the proposed method and suggesting some possible improvements, too.

METHODOLOGY
As typical of inverse gravimetric problems, the pro− posed algorithm can be divided into two main steps, the former consists in the data reduction until the isolation of the gravitational signal due to the particular mass anomaly to be investigated, the latter deals with the in− version of the isolated/residual field and the retrieval of some characteristics of the mass anomaly under study [Li and Oldenburg, 1998].

DATA REDUCTION AND SUBDUCTION EFFECT ISO-LATION
We will consider as observations a functional of the gravitational potential measured on a set of sparse points at satellite altitude and distributed all around the world (to fix the idea one can think to the second ra− dial derivative of the gravitational potential observed by the GOCE satellite). In order to fulfil the first step of the proposed procedure, i.e. to isolate the gravitational signal due to the investigated subduction zone, the gravitational effects of all the mass anomalies, apart from those of the subducting plate itself, should be modelled, preferably from gravity independent infor− mation, and then removed from the observations [Fors− berg et al., 2002;Sampietro et al., 2016; [Dziewonski and Anderson, 1981] but the density value of the oceanic crust from Carlson and Raskin (1984); units [g cm −3 ]. et al., 2018−a; Zaki et al., 2018−b]. In principle, this step leaves as unknowns only a set of parameters modelling the mass distribution related to the presence of the sub− ducting plate, thus reducing the space of the possible solutions, that is in general of infinite dimensions in the case of the inverse gravimetric problems [Sampi− etro and Sansò, 2012]. As said before, since the aim of this paper is to find a simplified way to model subduc− tion plates (from the gravitational point of view), we will just consider the crust−mantle density contrast ne− glecting any other existing crustal inhomogeneities. As well−known, this is a thorny issue due to the fact that such a kind of information is actually not available (at least on a global scale) with sufficient accuracy. Nonetheless, we start by modelling and removing from the observations the most external, and therefore prob− ably the most reliable, layers: namely, the gravitational effects of topography, bathymetry, and sediments. For this purpose, the ETOPO1 global relief model [Amante and Eakins, 2009] distributed by the National Geophysical Data Center (NGDC), with a spatial resolu− tion of 0.5°×0.5°, was used to compute the gravitational effects of land topography, ocean bathymetry, and ice surface (top of Antarctic and Greenland ice sheets) on a global spherical grid at a mean altitude of 10 km, i.e. on a so−called Brillouin sphere [Sansò and Sideris, 2013]. This altitude has been chosen since, in this way, the corresponding sphere is just outside the Earth's masses but it is close to them thus avoiding a too strong smoothing effect. The global gravitational effects have been computed by means of point−mass approximation as described in Reguzzoni et al. [2013], in details more than 180 million point masses have been used to pro− duce a global grid of the second radial derivatives of the gravitational effects of the ETOPO1 model at 10 km altitude with an expected accuracy better than 0.5 mE (1 mE = 10 −12 s −2 ). The same procedure has been applied to the CRUST1.0 sediment models [Laske et al., 2013]. Once computed, both signals have been removed from the actual dataset thus obtaining a reduced global grid of the second radial derivatives.
From this grid, a set of spherical harmonic coeffi− cients up to degree and order 360 were computed by performing a spherical harmonic analysis according to Colombo [1981] and Reguzzoni [2004].
At this point, we have removed the contributions of these layers from the GO_CONS_GCF_2_TIM_R5 GOCE global model [Brockmann et al., 2014] and we synthe− sized a grid of the second radial derivatives of the anomalous gravitational potential at GOCE mean satel− lite altitude, h = 250 km, from the residual coefficients. In this way, we obtain the gravitational effects of an Earth model "without" topography, bathymetry, ice− sheets, and sediments, in terms of the most informative GOCE observables. The residual gravity signal is basi− cally the sum of the effects of the subduction and of the crust−mantle interface. In particular, the latter term is expected, apart from anomalous regions, to be a smooth signal. The same holds also for the major density vari− ations in the lithosphere. The main task, at this time, is to stochastically exploit some properties of the residual field, such as its spatial correlation, in order isolate the gravitational effects of the subducting crust, which, vice−versa, is a sharp variation of mass density (see Figure 2), producing a quite sharp variation in the grav− itational field. Moreover, exploiting the statistical prop− erties of the reduced data decreases the signal that could be mistakenly considered as noise, in the sense of coloured noise, or as a residual effect of the isolated signal, i.e. due to the subduction plate [Braitenberg, et al., 2016;Sampietro, et al., 2017].
For this reason and in order to disentangle the differ− ent contributions inside the reduced gravitational signal, the procedure, summarised in Figure 4, has been applied as thoroughly explained in the sequel: 1) a mask is created within the region where the sub− duction zone is expected to be; 2) the empirical semi−variogram of the reduced grav− itational field is estimated using the observation points located close to the subduction region but yet outside the mask; 3) a Kriging procedure is used to predict the signal in− side the mask. The Kriging estimate [Wackernagel, 2013] allows ob− taining the gravitational field of a "regular" crust, i.e. a signal where the effects of the subduction are not evident.
The procedure was empirically tested on the Ker− madec−Tonga subduction zone. Basically, starting from the reduced second radial derivative of the anomalous gravitational potential, , synthesized on a regular grid with a resolution of 0.5° at 250 km altitude (see Figure 3−a) and considering the information coming from the SLAB1.0 model, the mask shown in Figure 3− b has been empirically defined.
The transformation from the geodetic, i.e. ellipsoidal, coordinate system to the geocentric system can be achieved using the following equations: (2) where, given a reference ellipsoid, a and b are respectively the semi−major and semi−minor axes, is the (dimensionless) first eccentricity, and N is the East−West radius of curvature.
The empirical semi−variogram γ( 12 ) was fit with a Steps for data reduction and isolation of the Tonga subduction zone gravitational field: (A) GOCE signal at satellite altitude reduced for the effects of topography, bathymetry, ice sheets, and sediments; (B) the region used for the semi− variogram computation excluding the points located within the mask in black; (C) the "regular" crust signal computed by Kriging to fill in the mask in (B); and (D) the isolated Tonga signal computed by subtracting the signal in (C) from the reduced GOCE signal in (A); units [mE].
theoretical one, in particular, many different models (e.g. bilinear, circular, spherical, pentaspherical, expo− nential, Gaussian, Whittle, Stable, and Matern [Wack− ernagel, 2013]) were tested and the best fitting, namely the exponential one, was chosen (see Figure 5). Using this theoretical semi−variogram function, the gravita− tional signal has been predicted by means of a standard Kriging solution (3) where is the interpolated value at point P 0 , located in− side the mask, using the measured values , i.e. grav− ity second derivatives with respect to the direction, namely East, North, and radial components, at points , where = 1,2,3... , in which the number of points in− side the mask is much smaller than the number of ob− servations , and is the weight of the measured value at the th location. The weights are obtained by solving the system reported in (4): where is a Lagrange multiplier used for the estimation error minimization, γ( ) is the value of the semi−var− iogram model for the distance between the th and k th observations, and γ( ₀) is the value of the semi− variogram model for the distance ₀ between the th observation and the point ₀ inside the mask. For de− tails on the Kriging solution, the interest reader is re− ferred to Wackernagel [2014]. Finally, the effect of the subduction is isolated by removing this signal from the originally synthesized (Figure 3−d). Even though the topography/bathymetry gravitational effects were removed, a high−frequency signal in correspondence with the Tonga−Kermadec trench and arc is still well−visible, negatively affecting the subsequent inversion. This is probably because of the presence of non−regular structures of the crust in correspondence of subducting plates, such as the ac− cretionary wedge or the magmatic arc, that are not taken into account in the reduction procedure. For ex− ample, Crawford et al. [2003] showed crustal thickness variations from more than 20 km to roughly 5 km over short spatial scales, which are obviously not included into our "regular" crust.
Note that this procedure could be directly applied to the observed gravity gradients instead of the synthe− sized grid of second radial derivative. However, in this way, one has to face the problem that along−orbit data are at different altitudes and therefore the semi−vari− ogram cannot, in principle, depend on the spherical dis− tance only.
In order to assess the validity of the proposed pro− cedure, a numerical test has been executed: basically, we randomly moved the defined mask in oceanic re− gions and we repeated the Kriging procedure. The basic idea is to verify the capabilities of the Kriging predic− tor to estimate the gravitational signal due to what we have called "regular" crust. This test shows that where the mask does not coincide with a subduction zone, the Kriging solution is able to predict the actual gravita− tional field with an accuracy of about 10.3 mE, in terms of standard deviation (STD). On the contrary, when the mask is located over a subduction zone, the residual signal shows a much higher STD value, e.g. over the Tonga subduction zone, the STD is 345 mE.
The whole procedure used to isolate the gravitational signal of the subduction plate is summarised in Figure  4, where it can be seen how the final signal used within the inversion is, basically, the difference between the observations reduced for the effects of the topography, bathymetry, ice model and sediments and the Kriging estimate of the "regular" crust.

INVERSION
The procedure to isolate the subduction signal ex− plained above is based on the use of the GO_CONS_GCF_2_TIM_R5 GOCE global model. How− ever, the aim of this work is to investigate the infor− mation content of the GOCE gravity gradients observed along the satellite orbit. For this reason, in order to "move" the estimated gravity effects of the "regular" crust into the sparse points along the orbit, the Kriging procedure has been applied all over the sphere and the resulting global grid has been analyzed in terms of spherical harmonic coefficients, as already done while computing the contributions of the topography, ba− thymetry, etc. Attention has to be paid to the fact that the results of this global Kriging interpolation represent the actual "regular" crust effects only within the study area, as shown in Figure 6. Of course, for the study of other subduction zones, the spherical harmonic coeffi− cients of the "regular" crust would be different.
Before applying the inversion, the observed gravity gradients are reduced for the effects of the "regular" crust as well as the topography, bathymetry, and sedi− ments, where their effects were synthesized from the corresponding spherical harmonic coefficients. In this way, the gravitational signal due to the density varia− tion between the subducting plate and the mantle is isolated from the rest of the observed signal.
Hence, it is possible to apply the inversion algorithm. Here, the main challenge is to find a methodology to dominate the well−known instability of such kind of inversion problems. In this sense, to further reduce the space of the possible admissible solutions, a simplified model that depends on a small number of unknown pa− rameters was required. Actually, one of the simplest models to describe a subducting plate is to model it as a dipping prism, which depends only on a set of eight unknown parameters (see Figure 7)  is the radius of the local sphere defined as: Note that, the three axes of the obtained Cartesian reference frame are almost pointing to the East, North and up direction, respectively.
In regards to the dipping prism model, it should be underlined that this approximation has a major model− ling advantage that allows an exact forward modelling of its gravitational field: in fact, a closed−form solution for the direct computation of the gravitational signal of a dipping prism was derived by Hjelt [1974]. However, in this work, we compute the gravitational effect using the point−mass approximation and we use the closed−for− mula only to check the accuracy of the forward model− ling. Later on, in fact, the point−mass approximation could allow modelling the gravitational effects of anom− alies depicted also with more complex shapes than a dipping prism, e.g. a dipping angle increasing with depth. Whatever forward modelling technique is used, the relationship between the unknown parameters and the observed gravitational field is strongly non−linear. Therefore, we decided to apply a SA procedure [Bertsi− mas and Tsitsiklis, 1993;Rossi et al., 2015] to find the set of parameters that best fits the isolated subduction signal (target signal). The Matlab implementation of the algorithm was used [Ingber, 1996]. Necessarily, the sys− tem must be initialized by imposing equal to a certain set ₀. Afterwards, a randomly chosen neighbourhood of ₀ (for instance, call it ₁) is drawn, the corresponding gravitational signal ( ₁) is computed and finally, the goodness of the solution is evaluated according to the following target function : where is the vector of the observed (target) second radial derivative of the reduced gravitational signal, ( ₁) is the predicted signal obtainable from the for− ward modelling of state ₁ of the SA algorithm, and is the observation error covariance matrix. According to the Metropolis−Hastings algorithm [Hastings, 1970], the state ₁ replaces the current state ₀ if the new state has a smaller value than the older one or if it passes an acceptance function. After that, the whole procedure is iterated until convergence is met.
The acceptance function and the size of the neigh− bourhood space, from which the new state is drawn, are in proportional to a given parameter (usually called "temperature"), which is decreased at each iteration ac− cording to a given cooling rule. The selection of the cooling function is an essential milestone in the per− formance of the SA algorithm. The decreasing temper− ature tends to force the current state towards the minima, moving only downhill. Decreasing the temper− ature too quickly could result in the state getting trapped in a local minimum, while decreasing it too slowly seems to waste the computational power through a very slow convergence [Kannan and Lakshmikan− tham, 2002]. The choice of the most suitable cooling rule for the above procedure is based on a closed−loop test that will be discussed in the next section.

CLOSED-LOOP TEST ON SYNTHETIC DATA
Before applying the proposed method to real GOCE data of a particular subducting plate, it was worthwhile to assess the algorithm and test its performance on a synthetic−realistic dataset. The closed−loop test basi− cally consists of building a synthetic subducting plate, with perfectly known geometry, simulating its gravita− tional signal, adding some noise, and applying the pro− posed inversion procedure. First of all, this permits to numerically assess the inversion strategy. Secondly, since the GOCE mission delivered six different kinds of along−orbit data, namely all the components of the Marussi tensor, the test could be also used to understand if some gravity gradients are more informative than oth− ers and, in case, to properly set the combination weights of all the different tensor components. In particular, we consider the three diagonal components of the Marussi tensor, namely , , and , adding to the observa− tions a 3 mE white noise, where , and stay for the East, North, and radial directions of the local triad. Again, we chose the Tonga subduction and the pa− rameters used to simulate the plate were picked from the literature [Nothard et al., 1996;Astiz et al., 1988] and reported in the first row of Table 1. The experiment has been repeated more than 100 times changing the realization of the noise in order to be able to statistically assess the performance. As for the cooling function re− quired by the SA procedure, here the Boltzmann cool− ing function [Geman and Geman, 1993] has been adopted, as follows: where ₀ is the initial temperature constant for = 0.
Among the different cooling rules, in fact, the Boltz− mann function led the system to the final state with the lowest value of the target function, consequently show− ing the best capability to converge towards the original parameters used to generate the simulated signal of the closed−loop test. Results in terms of STD of the differ− ences between the true parameters and the estimated ones are also reported in Table 1. Generally speaking, it can be seen that each gradient has different sensitivi− ties to the different parameters. In details, is highly− sensitive to all the parameters. On the other hand, regarding the and components, both show higher sensitivities to the parameters related to the per− pendicular direction. For example, gives a better prediction for the length in the X direction (Δx) and the strike angle ( ) than , which, in turn, is more sensi− tive to the width in the Y direction (Δy), the dipping angle ( ), and the y coordinate of the top face center than . Obviously, the aforementioned sensitivities highly−depend on the location and the geometry of the investigated subduction zone, for instance, for the strik− ing angle ( ), it is evident that a rotation of 90° of the dipping prism around the Z axis will basically exchange the rules of and .
Starting from the above results, it was possible to build a combined solution able to optimally exploit the information coming from the three different observ− ables of the main diagonal. In this research, the final solution is computed as a weighted average of the three individual solutions with the weights obtained for each unknown parameter defined as: where ( ) is the STD of the results obtained using the gravity gradients with respect to the direction for the th parameter of the dipping prism according to the ordering used in Table 1.

THE KERMADEC-TONGA SUBDUCTION CASE
The Kermadec−Tonga region extends from about 170° E to 170° W and 40° S to 10° S. It represents an ideal natural laboratory to study the presented inversion scheme for several reasons. The existence of a subduc− tion zone has been thoroughly studied using various seismological investigations [e.g., Oliver and Isacks, 1967;Mitronovas et al., 1969;Sykes et al., 1969;Mitronovas and Isacks, 1971;Barazangi et al., 1972]. The geometry of this subduction is fairly simple and well−studied [Bowman and Ando, 1987]. A large por− tion of the world deep earthquakes occurs in the sub− ducted slab in the Fiji deep seismic zone, and islands are located close to the deep earthquake epicentres, which allows recording in the shear−wave window.
From the gravitational point of view, this region is characterized by the presence of a series of geophysical signals at very different frequencies [Billen et al., 2003;Conder and Wiens, 2006], where the very long−wave− lengths (of the order of a few thousand kilometres) are dominated by the density heterogeneities inside the mantle. Superimposed on this long−wavelength signal, the Kermadec−Tonga and New Hebrides Trenches are marked by a high frequency (short−wavelength) signal. The observed bathymetry is dominated by short−to−in− termediate−wavelength (500 to 1000 km) signals caused by the crustal thickness variations, including the con− ( )= After checking the performance of the implemented algorithm by the closed−loop test, the synthetic signals have been replaced by the real ESA−GOCE along−orbit gravity gradients dataset and a solution for the Tonga subduction has been computed. In particular, the EGG_TRF_2 dataset [Bouman, et al., 2011] has been used instead of the original EGG_NOM_2 data [Rummel, et al., 2011]. For these data, a white noise has been as− sumed in the minimization of the three target functions (see (7) for the component). The inversion has been done targeting three individual solutions, namely a so− lution for each of the diagonal component of the Marussi tensor. Once computed, these three individual solutions have been merged by means of a weighted average using the weights resulted from the closed−loop tests (see (9). This procedure, namely the computation of a solution for each diagonal component of the Marussi tensor and the subsequent combination, allows to in− crease the observation number and to execute in par− allel the three SA optimizations.
The signal due to the subduction, as isolated by the Kriging procedure, is shown in Figure 3−d. It can be seen how the proposed algorithm allows to basically re− move the middle−low frequencies dominated by crustal thickness variations and mantle heterogeneities, thus highlighting the effects of the subduction (it is impor− tant to underline that the gravitational effects of to− pography, bathymetry and sediment layers were removed before applying the Kriging procedure). As al− ready stated before, a high−frequency signal that cor− responds to the Tonga−Kermadec trench and arc is not fully removed and it will be misinterpreted as the sig− nal of the subduction plate. However, it should be ob− served that this high−frequency signal is expected to have just minor effects on the estimation of the param− eters of our simple dipping prism model. The parame− ters estimated by inverting the three diagonal components of the Marussi tensor, namely the , , and , as well as the combined solutions, are reported in Table 2.
It can be noted that there is a good agreement, in general, among the three different solutions for Δx, Δy, Δz, x, and y with differences smaller than 10% with re− spect to the combined solution. As for the other pa− rameters, namely Δ and , it seems that the solution significantly deviates from the other ones with differences going up to about 60% with respect to the density contrast, 40% for the strike angle , and 15% for the dipping angle . This can be ascribed to the geometry of the subducting plate that mainly develops in the North−South direction, making the signal smaller than the other components. As a consequence, the impact of the non−perfect isolation of the subduct− ing plate as well as the simplification introduced by modelling the plate with a single dipping prism are am− plified, thus degrading the quality of the inversion of the component.
In any case, these results also show how the use of the three different functionals of the gravitational field makes the final solution more robust since the effect of possible mis−modelling acts differently on the compo− nents of the Marussi tensor. The use of these thre sec− ond derivatives of the gravitational potential produces a final result that is in a good agreement with the lit− erature values. For instance, the depth value Δz of the Tonga plate ranges between 670 km in Nothard et al. [1996] and about 700 km in the SLAB1.0 model [Hayes et al., 2012], which is in a good agreement with the 702 km of the combined solution. Similarly, the dipping angle , which is 40.6°, is in a good agreement with the values reported by Astiz et al. [1988] Nothard et al. [1996], and (2) according to Astiz et al. [1988]. model (which range between 40° and 55°). Also, an agreement with the SLAB1.0 model is found for the other planimetric parameters (namely x, y and ). As for the remaining parameters, i.e. the density contrast Δ and the length of the subducting plate Δx, since they cannot be observed by seismic data, no values have been found in the literature for a possible comparison.

DISCUSSION AND CONCLUSIONS
Although the complex shape of a subduction zone is highly simplified by a dipping prism or by a set of point−masses filling this dipping prism, the main esti− mated shape parameters, i.e. the planar extensions in the X and Y directions, prism depth, and both the dip− ping and strike angles, are in a relatively good agree− ment with the literature values. The overall differences between the obtained solutions and the seismic derived models with respect to the main planar parameters are less than 10%. The use of all the Marussi tensor diago− nal components, in principle, gives a more robust solu− tion due to the different sensitivity of each component towards each single model parameter.
Moreover, the results show how, with the advent of GOCE observations, gravity can give new useful infor− mation to improve our knowledge about subduction plates and tectonic dynamics, thus encouraging this kind of activities. In fact, regardless the ill−posed char− acteristics of the gravitational field inversion, an ade− quate gravity−related inversion technique would allow to retrieve information on the plate's density contrast with respect to the upper mantle, which cannot be re− covered from seismic data.
The SA is capable to converge and estimate the set of parameters generating the target signal with an accept− able error margin. For future research, the gravity in− version algorithm could be improved thanks to the SA capabilities of including numerical and physical con− straints and implementing hybrid functions to improve the search of the global minima of non−linear functions.
The information content of the different GOCE along−orbit components is similar. Consequently, the use of a target function that jointly minimizes the three main components of the Marussi tensor, instead of a posteriori combining them, could lead to the estimation of a slightly different set of parameters.
This paper is largely a proof−of−concept for using the GOCE data to study subducting plate structures. Naturally, the inversion of a simple uniform density prism adds very little to the understanding of the Tonga subduction plate and does not fully exploit the high quality of the GOCE data. A possible future perspective could be to break the prism into several sub−prisms along the strike, thus increasing the number of param− eters to solve, but offering the possibility of improving our understanding of the Tonga subduction complex and other plates more than a simple block prism ever could.