The static and time-dependent signature of ocean–continent and ocean–ocean subduction: the case studies of Sumatra and Mariana complexes


 The anomalous density structure at subduction zones, both in the wedge and in the upper mantle, is analysed to shed light on the processes that are responsible for the characteristic gravity fingerprints of two types of subduction: ocean–continent and ocean–ocean. Our modelling is then performed within the frame of the EIGEN-6C4 gravitational disturbance pattern of two subductions representative of the above two types, the Sumatra and Mariana complexes, finally enabling the different characteristics of the two patterns to be observed and understood on a physical basis, including some small-scale details. A 2-D viscous modelling perpendicular to the trench accounts for the effects on the gravity pattern caused by a wide range of parameters in terms of convergence velocity, subduction dip angle and lateral variability of the crustal thickness of the overriding plate, as well as compositional differentiation, phase changes and hydration of the mantle. Plate coupling, modelled within a new scheme where the relative velocity at the plate contact results self-consistently from the thermomechanical evolution of the system, is shown to have an important impact on the gravity signature. Beyond the already understood general bipolar fingerprint of subduction, perpendicular to the trench, we obtain the density and gravity signatures of the processes occurring within the wedge and mantle that are responsible for the two different gravity patterns. To be compliant with the geodetic EIGEN-6C4 gravitational disturbance and to compare our predictions with the gravity at Sumatra and Mariana, we define a model normal Earth. Although the peak-to-peak gravitational disturbance is comparable for the two types of subductions, approximately 250 mGal, from both observations and modelling, encompassing the highest positive maximum on the overriding plates and the negative minimum on the trench, the trough is wider for the ocean–ocean subduction: approximately 300 km compared to approximately 180 km for the ocean–continent subduction. Furthermore, the gravitational disturbance pattern is more symmetric for the ocean–ocean subduction compared to the ocean–continent subduction in terms of the amplitudes of the two positive maxima over the overriding and subducting plates. Their difference is, for the ocean–ocean type, approximately one half of the ocean–continent one. These different characteristics of the two types of subductions are exploited herein in terms of the different crustal thicknesses of the overriding plate and of the different dynamics in the wedge and in the mantle for the two types of subduction, in close agreement with the gravity data.


M O D E L S E T U P
The dynamics of the crust and mantle system are governed by the continuity, momentum and energy equations, which can be expressed as follows: where u is the velocity, p is the pressure, τ is the deviatoric stress, ρ is the density, g is the gravity acceleration, c p is the heat capacity at constant pressure, T is the temperature, K is the thermal conductivity and H is the total internal heating per mass unit.
Density varies with composition, temperature and pressure. To account for compressibility, we use the extended Boussinesq approximation (Ismail-Zadeh & Tackley 2010, and reference therein), and its implementation is described in detail in Appendix A. Here, we summarize only the main characteristics of this approximation. The balance equations are numerically integrated via the 2-D finite element (FE) thermomechanical code SubMar (Marotta et al. 2006) in a 2-D 1000 km wide and 700 km deep rectangular domain (Fig. 2a). During the progression of subduction, erosion/sedimentation, mantle hydration/dehydration and phase changes are allowed and implemented as detailed in Appendices A-C. The domain is discretized by a non-deforming irregular grid composed of 6813 quadratic triangular elements and 13 846 nodes, carrying a denser nodal distribution near the contact region between the plates, where the most significant gradients in temperature, pressure and velocity fields are expected. The sizes of the elements vary from 30 to 2.5 km, and the smallest elements are located in proximity to the active margin region to a depth of 300 km.
The marker-in-cell technique has been used to compositionally differentiate different types of materials (e.g. sediment, water, crust and mantle), and at the beginning of the evolution, 455 489 markers are spatially distributed with a density of 1 marker per 0.25 km 2 , identifying through a specific index the material particles of air, water, upper and lower oceanic crust, upper and lower continental crust and the mantle belonging to the upper plate. During the evolution of the system, each marker is advected using a 4th-order (both in time and in space) Runge-Kutta scheme. Based on the number and type of markers within each element, at each time it is possible to define the elemental composition and, consequently, all the properties, as specified in eqs (A2)-(A4) of Appendix A. The initial thickness of the oceanic crust is fixed at 10 km, subdivided into a 5-km-thick upper crust and a 5-km-thick lower crust, while the continental crust is 30 km, subdivided into a 20-km-thick upper and a 10-km-thick lower crust. To allow for topographic variations a 7.5-km-thick layer of sticky air has been introduced above the crust. Only in the oceanic part of the domain, a sticky water layer of thickness h w is sandwiched between the sticky air and the crust (Fig. 2a). When we refer to topography, we do not mean the upper boundary of the model, which does not deform, but the surface defined through the envelope of the topmost continental and/or oceanic crustal markers. Several comparative studies have demonstrated that, when appropriate low density and low viscosity values are considered (as low as 0-1 kg m −3 for density and 10 18 −10 20 Pa·s for viscosity), by means of the sticky air method the air/crust interface behaves in a manner that is very similar to that of a true free surface (Schmeling et al. 2008;Crameri et al. 2012;Ruh et al. 2013). In this study we assume a value of 10 18 Pa·s for the air/water viscosity and values of 1 and 1000 kg m −3 for the air and water density, respectively.
The model combines a linear viscous rheology for the sublithospheric mantle with a linear viscoplastic rheology for the lithosphere, as detailed in Appendix A. The initial global thermal configuration of the system corresponds to a conductive gradient throughout the lithosphere, with a temperature that ranges from 300 K at the surface to 1600 K at its base and a uniform temperature of 1600 K below the lithosphere to the bottom of the model. At the beginning of subduction, the base of the thermal lithosphere is fixed at a depth of 80 km under the lower plate and at 80 or 150 km under the upper plate (red solid and dashed lines in panel a of Fig. 2). A thickness of 80 km corresponds to either an oceanic  (6) and (7). Distances are not to scale.
lithosphere of approximately 40 Myr (Turcotte & Schubert 2002) or to a thinned continental passive margin extending at a medium-to-slow spreading rate of 2-3 cm yr -1 (e.g. Marotta et al. 2018). The 1600 K isotherm defines the base of the thermal lithosphere throughout the evolution of the system.
We use different type of boundary conditions, either Dirichlet or Neumann type, along different boundary of the 2-D domain (Fig. 2a), where we prescribe values of velocities and temperatures or their normal derivative.
The thermal boundary conditions correspond to fixed temperatures at the upper (300 K) and lower (1600 K) boundaries of the model and to zero thermal flux through the vertical sidewalls. To minimize the potential effects of an unrealistic motion of the two plates with respect to the deep mantle, our model setup is constructed such that the trench is fixed with respect to the bottom of the 2-D convective cell. This assumption, which is also encountered in the perspective of the successive application to the case studies of the Sumatra and Mariana subduction complexes, is consistent with the absolute plate motion by Wang et al. (2018) at these two locations. Table 1. List of the models implemented for the present analysis. U cc /L cc : stratified (yes) or homogeneous (no) continental crust of the upper plate; u s : prescribed subduction velocity; θ s : prescribed subduction dip angle; c f : degree of plate coupling; L c : initial thickness of the upper plate; h t : maximum topographic height before the activation of erosion; c fo : degree of coupling between the subducting plate and the sticky air layer. ms: modified sedimentation algorithm. Moreover, no-slip conditions are assumed along the upper (top of the air/water layer) and lower boundaries of the 2-D domain and along the right vertical side, with the exception of the uppermost 80 km, where zero normal viscous stress is assumed. Zero normal viscous stress is also assumed along the vertical left-hand side (Fig. 2a). To force the subduction of the oceanic lithosphere, a convergence velocity u s is prescribed along the bottom of the oceanic crust. Velocities of 2, 4, 5 and 8 cm yr -1 have been chosen, which are representative of slow, intermediate and fast subductions, respectively. The same velocity is also imposed along a plane that extends from the base of the oceanic lower crust, behind the trench, to a depth of 80 km, with an angle θ s of 30 • , 45 • and 60 • (dashed black arrows in panel a of Fig. 2). These dip values allow the development of low-, medium-, high-dip subductions and are based on the compilation by Syracuse & Abers (2006) who averaged the slab dips between 50 and 250 km depth, for all the known subductions. Flow is allowed across the left edge to guarantee that the total mass is conserved in the cell.
Finally, a different plate coupling along the contact plane between the lower and upper plates is implemented, modifying the classical split-node technique (Jungels 1973;Jungels & Frazier 1973) by introducing a plate coupling factor c f that varies from 0 to 1 to indicate a shear-free or a fully coupled contact surface, respectively (Fig. 2b). The details about the implementation of the plate coupling factor c f are presented in Appendix D.
The assumption of a total degree of coupling (c f = 1) corresponds to a condition of continuity in the velocity field through the contact surface between the two interacting plates, that is commonly used in similar studies, even when a low viscosity channel or a weak zone are implemented to ensure the separation between the two plates (e.g. Billen & Hirth 2007;Duretz et al. 2011) intrinsically contains the condition of continuity. In our study, instead, we simulate a varying degree of coupling in terms of a discrete structural discontinuity. We assume a reference model that simulates a fully coupled subduction in order to enlighten its effects on the gravitational signature of subduction. Table 1 lists the models implemented in this study and their main characteristics.
Because we are interested in estimating the perturbation induced on the gravity field only by subduction, for each element e of the numerical grid, at each time we calculate the density anomaly as follows: where ρ e (t) is the density within element e at time t and ρ e (t 0 ) is the density within the same element at the beginning of subduction. The elemental density anomaly ρ e (t) computed by eq. (4) are then used to calculate the corresponding gravitational contributions at the observational point S, g(S, t), as follows: where G is the universal gravitational constant, 6.67 × 10 −11 m 3 kg −1 s −2 ; b is the distance between the observational point S and the centre of the volume element d e inside grid element e; e is the volume of the prism, infinitely extended along the direction perpendicular to the model, whose section coincides with the grid element (Fig. 2c); nelem is the total number of the grid elements. After few mathematical steps eq. (6) becomes where A e is the area of the grid element, x b e and y b e are the coordinate of its barycenter and x S and y S are the coordinates of the observational point S.

M O D E L R E S U LT S A N D D I S C U S S I O N
The description of the results concerns the general thermal state of the system and the dynamics in the wedge area, the density distribution and the gravity changes induced by subduction. First, we illustrate the results of a reference model (OC 2 model in Table 1) (Section 3.1). This model simulates a fully coupled ocean-continent subduction (c f = 1), with a subduction velocity u s = 5 cm yr -1 , a subduction dip angle θ S = 45 • . Secondly, we compare the results of models characterized by different parameters in order to enlighten the effects of a variation of the degree of coupling c f (1, 0.5, 0.25), the subduction dip angle θ S (45 • and 60 • ), the subduction velocity u s (2, 4, 5 and 8 cm yr -1 ) and the environment (ocean-ocean and ocean-continent) (Section 3.2).

General thermomechanics
The general thermomechanical behavior of the crust-mantle system during subduction has already been extensively described in the literature (e.g. Marotta et al. 2007;Billen 2008;Afonso & Zlotnik 2011;Roda et al. 2012;Gerya 2015;Regorda et al. 2017;Dai et al. 2018;Stern & Gerya 2018;Wang et al. 2019). Here, we will focus on those aspects that we consider to be related to the gravitational signature and its changes in time, which are the objectives of this work. Fig. 3 shows the temperature and velocity fields predicted by the reference model OC 2 at different times after the beginning of subduction. The system is characterized by the following: (i) The sinking of the cold lithosphere into the hot mantle and its progressive steepening at great depths (Fig. 3d); (ii) A convective flow that is activated above the lower plate during the early stages of subduction (Figs 3a and b) and progressively enlarges with the sinking of the slab (Figs 3c and d); (iii) The large-scale flow that favours the rising of hot mantle in the wedge area, producing a progressive, but local, thermal thinning of the overriding plate (Figs 3c and d).
Focusing on the mass redistribution shown in Fig. 4, it is also worth highlighting the following: (i) The subducted crustal material originates either from the lower plate or from the upper plate, from which it is scraped (Fig. 4a); (ii) The scraping effect is maximum within the first 5 Ma (Figs 4a and b); (iii) During the early stages of subduction, crustal erosion and thinning are controlled by the mechanical coupling between the interacting plates (Figs 4a and b); (iv) In the advanced stages of subduction, however, the main role of crustal erosion and thinning is played by the local mantle convective flow in the wedge area (Fig. 4c); (v) Since the beginning of subduction, a widening of the trench trough occurs until a stationary configuration is reached after a few million years (cyan colour in   (vi) After a few million years, the erosion/sedimentation mechanism is activated, with a vast amount of sediments that, with water, contributes to filling the trench and that sinks to great depths (black dots in Figs 4a-c); (vii) The p-T conditions become favourable for activating mantle hydration and serpentinization in the wedge area within the first 2 Myr from the beginning of subduction (yellow area in Figs 4a-c).

Density anomalies
The density anomalies predicted by the OC 2 model at different times after the initiation of subduction are shown in the insets of Fig. 4.

796
A.M. Marotta et al. Alternating regions of negative (blue colour) and positive (yellow to red colours) density anomalies are localized at shallow depths around the arc-trench region during the whole evolution (insets in Figs 4a-c); a dominant positive density anomaly, increasing in time, characterizes the core of the subducting lithosphere to great depths (Figs 4d and e). Both the extension and the magnitude of the density contrasts are sensitive to the age of subduction.
The negative density anomalies can be attributed to the following: (i) The negative topography associated with the oceanic trench. The depression that constitutes the oceanic trench is filled with water and sediments and produces the most intense (and shallowest) negative anomalies (insets of Figs 4a-c); (ii) The light material entrapped in the subduction complex, particularly sediments and slices of continental crust scraped by the upper plate. The sediments and the continental material play a peculiar role: if their presence is initially important (Fig. 4a), then at longer times, their concentration decreases because they are transported at depth by the flow generated by the descent of the subducted plate (Figs 4b and c); (iii) The hydration of the mantle wedge leading to the development of serpentine (yellow colour in Figs 4a-c) that has a density that is approximately 200 kg m -3 lower than that of the non-hydrated mantle; (iv) The general warming of the system at the mantle levels of the wedge area (see the deflection of the isotherms 800 and 1600 K in Figs 4b and c).
The positive density anomalies can in turn be attributed to the following: (i) The sinking of cold and dense lithospheric mantle into warm and light material; (ii) The phase transitions within the subducting plate triggered by the increase in pressure p and temperature T. In fact, mineralogical phases that are stable at lower depths are replaced by denser phases that are stable at higher p-T conditions. The effect of phase transitions adds to the thermal contraction; (iii) The mantle that rises and fills the space previously occupied by the crust that has been scraped from the overriding plate (Figs 4a-c); (iv) The general cooling and thermal contraction of the system far away from the trench.

Gravitational contribution
The profiles of the gravitational contributions, calculated for the OC 2 model using eqs (5) and (6), at time intervals of 1 Myr since the beginning of subduction and along the top boundary of the model domain are shown in panel a of Fig. 5. Each profile is characterized by two maxima and a minimum within a distance of 200 km (Fig. 5a). The most intense maximum is located on the upper plate, while the secondary one is located behind the trench. The local minimum is positioned approximatively above the trench. The gravity anomalies decrease progressively outwards, assuming an almost constant value in the far field. With the progression of the subduction, the intensity of the gravity anomalies increases globally, accompanied by a coeval increase in the difference between the absolute maximum and the minimum. For the OC 2 model, the positions of the two maxima and of the minimum remain almost unchanged for the entire evolution. The analysis of the separate contributions of the positive and negative density anomalies (Figs 5b and c, respectively) allows us to clarify the origin of the peculiar style of the gravitational contribution and its variation in time and space. The bulk of the gravity change due to the negative anomalies in the surrounding of the trench occurs within the first two million years ( Fig. 5c) and is attributable to the formation of the trench trough, to the subduction of crust and sediments and to the mantle wedge hydration. Subsequently, the most significant changes occur far from the subduction complex, and these changes can be associated with the progressive but widespread cooling of the system. Conversely, the progressive growth of the positive density anomalies (Fig. 4) induces an equally continuous increase in their gravitational contribution, although with different rates (Fig. 5b). The following four phases in the time variation of the gravity contribution can be distinguished: (i) A first phase, from the beginning of subduction to approximately 5 Myr. This phase is characterized by the highest rate of growth in the surroundings of the trench, with a total variation of approximately 200 mGal (blue lines in Fig. 5a). During this first period, the subducted slab reaches a depth of more than 200 km (Fig. 3b), and most of the negative and positive density anomalies are localized at shallow/intermediate depths around the trench (inset of Fig. 4b); (ii) A second phase, from 5 to 9 Myr (cyan to yellow lines in Fig. 5a). During this phase, the gravitational contribution continues to increase, but at a lower rate, and the positive density anomalies still grow but at great depths (Fig. 3c, and inset of Fig. 4c); (iii) A third phase, from 9 to approximately 11 Myr (yellow to dark orange lines in Fig. 5, panel a). This phase is characterized by a gravity change of approximately 60 mGal within just 2 Myr. Note that this effect occurs when the subducted slab crosses the olivine − spinel transition at a depth of 410 km (Fig. 4, panels d and e), where a density change of approximately 100 kg m -3 occurs. The great depth at which the density variation occurs makes the corresponding perturbation induced on the gravity field have a regional character; (iv) A fourth phase, from 11 Myr onwards (dark orange to pink lines in Fig. 5a), during which the growth rate is stabilized everywhere at low values of the far field. Fig. 6 shows the rate of change of the gravitational contribution of the density anomalies predicted by the OC 2 model. The plotted rate of change corresponds to the average value calculated over a time span of 0.5 Myr, which confirms that during the first phase, the gravitational contribution of the negative density anomalies increases rapidly at a rate as high as 0.05 μGal yr -1 (Fig. 6c). Afterwards, the rate of change remains smaller than 0.01 μGal yr -1 . Regarding the gravitational contribution of the positive density anomalies (Fig. 6b) during the first phase, the corresponding rate of change is comparable to that associated with the negative density anomalies and is reduced by almost a factor of two during the second phase. The discontinuity observed at approximately 9 Myr from the beginning of subduction (yellow lines in Figs 6a and b) is evident here in terms of a rate of change varying from approximately 0.02-0.03 μGal yr -1 , moving from the lower to the upper plates; values close to 0.04 μGal yr -1 characterize the surroundings of the trench (Figs 6a and b). The low values of the rate of change characterizing the mature phases indicate that, at long times, the contribution from the positive density anomalies becomes quasi-stationary.

Coupling factor c f
To investigate the effects of a different degree of plate coupling on the gravitational field, we implemented a different set of models with the same values for all the characteristic parameters, with the exception of the value of the coupling factor c f , set equal to 1 (total coupled system), 0.5 and 0.25 (almost completely decoupled system). Fig. 7 shows how a different plate coupling degree along the contact between the subducting and overriding plates impacts the global thermomechanics of the subduction complex, with panels a and b corresponding to the fully coupled OC 2 model and panels c and d corresponding to the OC 5 model, characterized by a c f = 0.25. During the early stages of subduction, no significant differences in the global thermal field are evident, with the exception of the slightly lower temperatures in the wedge area occurring as the plate coupling decreases (compare panel c to panel a of Fig. 7). Over time, the difference in the thermal state is intensified and, as the plate coupling decreases, the cooling of the wedge area progressively increases (compare panel d to panel b of Fig. 7) while, at great depths, the subducted slab remains warmer (compare panels e and f of Fig. 7, corresponding to c f = 0.25, to panels c and d of Fig. 3, corresponding to c f = 1). The thermomechanics remain unchanged at great distances from the trench.
The explanation of the enlightened changes in the thermal field stands on the more vigorous small-scale convection when the plate coupling is higher in the corner region between the subducting and overriding plate, which enhances the entering of cold material into the mantle due to the active subduction, making the fully coupled OC 2 model colder at depth compared to the decoupled OC 5 model. At the same time, as already stated in Section 3.1.1, the large-scale flow favours the rising of hot mantle in the wedge area, producing a progressive local thermal thinning of the overriding plate (Figs 3c and d). The concomitance of these two factors makes the system warmer at shallow depths and colder at deep depths. For low coupling (Figs 7c-f) the two plates are partially or totally mechanically disjoined and the low intensity convective cornel flow is not able to perturb the thermal field in the wedge area, where the thermal state becomes increasingly depressed over time. At the same time, the progressive cooling of the wedge area further inhibits the convective flow in the wedge which in turn disfavors the subduction of cold material for the lower plate coupling model, which thus remains warmer at great depths.
The role of plate coupling appears to be even more evident from the distribution of the markers in the wedge area ( Fig. 8). At the early stages of subduction, the scraping of the crustal material belonging to the upper plate is increasingly less effective as the plate coupling decreases (Figs 8a and b, compared to Fig. 4a). A smaller amount of the scraped crust is subducted at great depths over time compared to the OC 2 model (Fig. 4c), and over long times, the recycling flow of crustal material inside the hydrated mantle wedge is more intense (Figs 8c and d). Focusing on this last aspect, the maximum depth of the hydrated area is greater for strong plate coupling ( Fig. 4 with respect to Fig. 8); conversely, the horizontal extension of the hydrated area increases as the coupling decreases, favoured by the low temperatures that characterize the wedge area.
The differences in the density anomalies predicted by models characterized by different couplings (insets of Fig. 8) increase with the age of subduction. The main differences occur at low depths, where the less effective mantle wedging occurring for low plate coupling causes a shallow positive density anomaly of smaller intensity (insets of Fig. 8 compared to insets of Fig. 4). The positive density anomaly is instead mainly due to the ascent to low depths of the lower crust, which is heavier than the upper crust. The negative density anomaly, associated with the trench trough, the subducted crust and sediments and the mantle hydration, increases and widens laterally as the plate coupling decreases (Figs 8c and d,compared to Fig. 4c). Finally, the diffuse negative density anomaly associated with the warming of the mantle in the wedge area disappears for low values of plate coupling. Fig. 9 shows the gravitational contribution predicted by the OC 3 and OC 5 models, accounting for plate couplings of c f = 0.5 and c f = 0.25, respectively. The decrease in the maximum value of the positive contribution accompanying the decrease in plate coupling is due to both the warming of the internal part of the subducted plate, which is lighter for low plate coupling with respect to the strong coupling (Figs 7e and f), and to the missing wedging mantle flow (Figs 8b and d). In contrast, the gravitational contribution due to negative density anomalies increases both in magnitude and in width as the plate coupling decreases, reflecting the differences in the negative density anomaly distribution discussed above and particularly the progressive enlargement of the hydrated mantle area where a very intense recirculation of crustal material occurs (Figs 8b and d).
Note that with the progression of the subduction, the amplitude of the gravity trough increases as the plate coupling decreases. Although the degree of plate coupling does not affect the position of the local minimum that remains localized above the trench throughout the evolution, the position of the absolute maximum moves forward over the upper plate as the degree of plate coupling decreases ( Fig. 9), which is correlated with the main role played by the heavy subducted slab in localizing the bulk of the positive density anomaly further away from the trench, below the upper plate, as subduction progresses. Fig. 10 shows the temperature and velocity fields (panels a and b), the marker distribution (panels c and d) and the density anomalies (insets of panels c and d) predicted after 2.5 and 13 Myr from the beginning of subduction by the OC 7 model that is characterized by a subduction velocity u s = 5 cm yr -1 , a coupling factor c f = 0.5 and a subduction dip θ s = 60 • .

Subduction dip θ S
With respect to the homologous 45 • OC 3 model, a higher subduction dip angle allows the development of a more intense mantle wedge flow, responsible for a stronger erosion and thinning of the upper plate, with a larger amount of light material (both crust and sediments) dragged at great depths (compare panels c and d of Fig. 10 to panels a and c of Fig. 8). The local thermal warming occurring at mantle levels in the wedge area also generates unfavorable p-T conditions for mantle hydration, which in turn limits the recirculation of light material. All these factors generate intense and localized negative and positive density anomalies, as shown in the insets of panels (c) and (d) of Fig. 10. The substitution of the scraped crust and the rising of cold mantle is responsible for shallow positive density anomalies that are considerably more intense than that predicted by the homologous 45 • OC 3 model (compare panel d of Fig. 10 to panel c of Fig. 8). Conversely, at greater  Regarding the negative density anomalies, the component associated with the trench trough, the subducted crust and sediments and the mantle hydration remains almost unchanged in magnitude but assumes a thinner shape along the subduction plane. Furthermore, a diffuse negative anomaly appears at the mantle level of the wedge area, similar to that predicted by the fully coupled OC 2 model (Fig. 4c). However, it is now less intense because the lower plate coupling makes the local crustal erosion due to the mantle wedging flow not accompanied by coeval and equally effective local heating (compare isotherms 800 and 1500 K in Fig. 10d, and Fig. 4c). Fig. 11 shows the gravitational perturbation due to the density anomaly shown in Fig. 10. A decrease of approximately 20 mGal in the magnitude of the gravitational contribution due to the negative density anomalies occurs with respect to the 45 • OC 3 model during the whole evolution (compare panel c of Fig. 11 to panel c 1 of Fig. 9); an increase in the gravitational contribution due to the positive density anomalies also occurs, up to more than 50 mGal in the mature stages of subduction. In addition, during evolution, the lateral extension of the gravity trough decreases with time because the absolute maximum moves towards the trench as a consequence of the intensification and localization of the shallow positive density anomaly (insets of Figs 10c and d). This behaviour is opposite of that occurring with the 45 • OC 3 model. Finally, because the subducted plate reaches the olivine-spinel transition in a more vertical way, the jump in the gravitational contribution occurs earlier and more rapidly (compare Fig. 11b, to Fig. 9b 1 ).

Subduction velocity u s
To identify the effects induced by a different subduction velocity on the gravitational field, we compare the predictions from three models characterized by the average values of all the characteristic parameters, with the exception of the assumed subduction velocity, set equal to 8 cm yr -1 (OC 6 model), 5 cm yr -1 (OC 3 model) and 2 cm yr -1 (OC 14 model) (panels a i , b i and c i of Fig. 12, respectively). In this case, the comparison is not performed after the same time span since the beginning of the subduction but rather after a different time span for each model, when the same amount of ocean is consumed, specifically 200 and 600 km, corresponding to approximately 2.5 and 7.5 Myr for u s = 8 cm yr -1 , 4 and 12 Myr for u s = 5 cm yr -1 , and 10 and 32 Myr for u s = 2 cm yr -1 . The depth of 200 km is chosen because the subducted slab is well developed but still unaffected by the olivine-spinel phase boundary, while the depth of 600 km is chosen because it is comparable with the maximum length of the Benioff planes. The major effects of a change in subduction velocity are in the thermal field of the wedge area and of the subducted slab (Fig. 12). From the general thermomechanical perspective, the more efficient advection characterizing the high-velocity model OC 6 of panels a 1 and a 2 is responsible for a thinner and colder slab that favours maintaining the prescribed subduction dip angle at great depths (panel a 2 of Fig. 12, compared to b 2 and c 2 ). In addition, the thermal field remains almost unchanged far from the trench (Fig. 12). A decrease in the subduction velocity results in a thermal thickening of the subducted slab at shallow depths, in a more effective thermal erosion at high depths, in a significant variation of the deep dip and in a thermal thickening of the upper plate (Fig. 12, panels b 1 , b 2 , c 1 and c 2 ).
If we focus on the wedge area and on the distribution of the material particles ( Fig. 13, with the fast OC 6 model in panels a and c and the slow OC 14 model in panels b and d), it is evident that the thermal differences are not accompanied by equally significant differences in the distribution of material particles and by significant differences in the shallow density anomalies (compare panels c and d of Fig. 13 to panel c of Fig. 8).
Instead, at the deep lithosphere mantle level and at great depths, the thermal thickening characterizing the slowest subduction velocity model OC 14 induces a widespread positive density anomaly close to the lithosphere base and a positive anomaly associated with the subducted slab with a thickness of almost one and a half times (compare panels e and f of Fig. 13).
Fig. 14 shows the gravitational perturbation at times when the subducted slabs reach depths of 200 km (panels a 1 , b 1 and c 1 ) and 600 km (panels a 2 , b 2 and c 2 ) from the fastest to the slowest models, from top to bottom. Although the contribution of the negative density anomalies (dotted black lines) is not significantly influenced by the subduction velocity, the contribution from the positive density anomalies (dashed black lines) increases as the subduction velocity decreases due to the thermal thickening of the subducted plate and to the anomaly of the positive density that develops near the base of the lithosphere (also see the insets of panels e and f of Fig. 13). The latter, due to its widespread nature, also intensifies the asymmetry in the gravity profile, increasing the difference between the minimum value of the gravitational contribution and its far field value as the subduction velocity decreases. Finally, note that at any time, the width of the gravity trough increases as the subduction velocity decreases.
It is certainly interesting to compare the rate of change of the gravity anomalies for different subduction velocities. Based on the considerations made in the previous paragraphs regarding the variations in time of the areal distribution and of the magnitude of the density anomalies obtained for different subduction velocities, we expect that the rate of change of the gravity anomalies varies significantly as the subduction velocity varies. This is in fact what occurs in Fig. 15, where the rates of change for the OC 6 (u s = 8 cmyr -1 ) and OC 14 (u s = 2 cmyr -1 ) models are presented. As expected, the differences can be significant and the rate of change decreases to one third when the subduction velocity decreases from 8 to 2 cm yr -1 (Fig. 15), making slow subduction appear to be a quasi-static process even from the gravitational perspective. However, note that when the rates of change are normalised with respect to the corresponding subduction velocity, the graphs, with the exception of some differences due to local-scale features, become comparable (Fig. 16), making it possible to estimate a maximum rate of change equal to 0.008 μGal yr -1 for each cm/yr of subduction velocity.

Ocean-ocean context
Here, we illustrate only the major differences between predictions from ocean-ocean models and the corresponding ocean-continent models. Fig. 17 shows the marker distribution and the density anomalies (in the insets) predicted by the OO 1 , OO 4 and OO 7 models, with the coupling factor decreasing from top to bottom, from 1, 0.5 and 0.25, after 2.5 and 13 Myr from the beginning of subduction, at a velocity of 5 cm yr -1 , with a subduction dip angle of 45 • . With respect to the corresponding ocean-continent-type models (models OC 2 - Fig. 4, OC 3 and OC 5 - Fig. 8), a generally larger amount of light material (both crust and sediments) recirculates in the wedge where a wider hydrated area is further developed. A significant erosion of the crust belonging to the upper plate occurs, allowing the exhumation of the recirculated light material. This peculiar mass redistribution is at the origin of the very intense and shallow negative density anomalies that develop over an area wider than 60 km. Fig. 18 shows the gravitational contribution at 4 and 12 Myr from the initiation of subduction for the same models shown in Fig. 17, with the coupling factor decreasing from top to bottom, from 1, 0.5 and 0.25. Unlike the ocean-continent context, plate decoupling does not significantly affect the shape of the gravitational contribution pattern. The aspects that deserve to be noticed are a slight increase in the positive contribution accompanying the decrease in the coupling of the plates and a slightly more significant increase in the negative contribution with the increase in the coupling (compare panels a and c of Fig. 18). The most important characteristic is instead the greater width of the gravity trough, on the order of more than 200 km, that occurs in the ocean-ocean environments compared to 100 km in the ocean-continental environment (Fig. 14, panel b 2 ) due to the peculiar wider distribution of the recycled material (Fig. 17).

Tectonic setting
The Sumatran and Mariana subductions are considered to be two classical tectonic settings representative of an ocean-continent subduction and an ocean-ocean subduction.
The Indonesian Island of Sumatra represents an interesting natural laboratory for geodynamic and geophysical investigations because it is one of the most active regions in the world. In fact, numerous large earthquakes have occurred in the Sumatra subduction zone within the last two centuries, including the magnitude 9.3 Sumatra-Andaman earthquake on 26 December 2004 (Gamage 2017). The subduction zone of Sumatra results from the underthrusting of the Indo-Australian Plate subducting beneath the Sunda plate, Andaman and Burma microplate (Shapiro et al. 2008;Gamage 2017) along three main trenches, Java, Sunda and Andaman-Nicobar (Nielsen et al. 2004), with a convergence  The convergence is nearly orthogonal along the Java trench and subparallel along the Andaman Islands trench (Moeremans et al. 2014). The transition between the two subduction regimes occurs south of the Sunda Strait (Malod & Kemal 1996). Here, where the plate boundary geometry changes abruptly, tomography and seismic studies show slab tearing, with a distinct gap in seismicity between depths of 300 and 500 km beneath the central Sunda (Java) block (Kundu & Gahalaut 2011). In addition to strong strain partitioning, the area presents other peculiar tectonic settings, such as active spreading in the backarc beneath the Andaman Sea and lateral age variability of the incoming Indian Plate. Beneath Sumatra, a young oceanic lithosphere of approximately 40-60 Myr is subducting (Schluter et al. 2002;Shapiro et al. 2008), corresponding to a thickness of 80-100 km (Turcotte & Schubert 2002). In contrast, off Java and Andaman Islands, the age of the oceanic lithosphere varies from 70 to 120 Myr and from 70 to 90 Myr (Moeremans et al. 2014). This aspect influences the style of deformation and seismicity along the entire arc. Where the lithosphere is younger and consequently warmer and more buoyant, the subducting plate and the mantle wedge appear to be strongly coupled with a Wadati-Benioff zone of approximately 30-45 • (Shapiro et al. 2008;Kundu & Gahalaut 2011). In contrast, in adjacent areas to the north and south, the seismic zones have higher dip angles (Shapiro et al. 2008;Kundu & Gahalaut 2011). The variability in the convergence velocity, particularly in the component perpendicular to the trench, as well as in the dip angles, from 30 • to −45 • to higher angles, indicates that the appropriate velocity values of 8, 5 and 2 cm yr -1 (with 5 cm yr -1 as the most representative value) and the three dip angles of 30 • , 45 • and 60 • (with 45 • the most representative one) are appropriate for the modelling of the EIGEN-6C4 gravitational disturbance characterized by a 2-D dipolar pattern perpendicular to the trench, first requiring a 2-D modelling approach as performed herein to capture the characteristics of the gravity pattern.
Regarding the age of subduction, for these case studies, we take as a reference the youngest age of the Sumatran oceanic lithosphere of 40 Myr.
Taken collectively, within the frame of a 2-D approach used to model the 2-D pattern of the EIGEN-6C4 gravitational disturbance, similar features can be found in the Izu-Bonin-Mariana (IBM) trench. Stretching for 2800 km (Pearce et al. 2015) accommodates the oblique west-dipping subduction (Faccenda et al. 2018;Kong et al. 2018) of the Pacific beneath the Philippine Sea Plate at a velocity of approximately 9 cm yr -1 along the Izu-Bonin arc and 5 cm yr -1 along the Mariana arc (Kong et al. 2018). Since the Eocene, the plate-tectonic history of this region has been affected both by multiple clockwise rotation events of the Philippine Sea Plate, as evidenced by paleomagnetic data (Hall et al. 1995;Sdrolias et al. 2004), and by the northeast migration of the triple junction between the Eurasian, Philippine Sea, and Pacific Plates from 30 to 17 Myr (Castle & Creager 1999). Along with the triple junction motion, the northern Izu-Bonin trench retreated over 1000 km to the northeast, producing the opening of the Shikoku basin (between 25 and 15 Myr, Hall et al. 1995;Anderson et al. 2017 (van der Hilst & Seno 1993), corresponding to a thickness of 160 km (Turcotte and Schubert, 2002). Several studies of seismic tomography and seismicity have clearly defined the geometry of the subducting slab in depth. The variation in the geometry of the subducted slab along the arc is a peculiar aspect of the region. van der Hilst & Seno (1993) proposed that beneath the Izu-Bonin arc, the combination of young oceanic lithosphere and rapid trench migration produced a shallow subduction angle, with the slab deflecting horizontally at the 670 km discontinuity (van der Hilst & Seno 1993;Miller et al. 2005 In contrast, southward, below the Mariana trench, the slab plunges steeply into the lower mantle (Faccenda et al. 2018;Holt et al. 2018). The abrupt change from a horizontal to a vertical slab is accommodated by tearing at the junction of the Izu-Bonin and Mariana arcs at ∼26 • N, where the Ogasawara Plateau occurs at the junction (Anderson et al. 2017). In fact, tomographic and seismic images reveal a distinct change in the seismic property of the subducted Pacific plate beneath the Izu-Bonin arc, with a seismic gap at approximately 300-400 km, interpreted as a probable slab tear (Kong et al. 2018). A N-S trending slab tear is also suggested at the southern end of the Mariana arc, known as the Challenger Deep (Miller et al. 2006).
Due to the similarities in convergence velocity, dip angles and age, as well as their variability along the trench, we consider the Sumatra values for these parameters to also be appropriate for the Mariana case study, focusing on the effects of the most important difference between these two types of subduction, namely, the nature of the overriding plate, oceanic versus continental. Fig. 19 shows a regional map of the gravity disturbance in the surroundings of the Sumatra and Mariana subduction complexes based on the static global combined gravity field model EIGEN-6C4 (Förste et al. 2014) and evaluated in the spherical approximation at a height h S =5 km above the spheroid of radius R = 6378.137 km. We consider this altitude to smooth slightly the geodetic signal and make it more suitable to be compared with the modelled one. The geodetic functional gravity disturbance is defined as the difference between the magnitude of the gravity g and the magnitude of the normal gravity γ : OO 7 c f = 0.25 Figure 18. Solid lines indicate the gravitational contribution of the density anomalies predicted by OO 1 (panels a 1 and a 2 ), OO 4 (panels b 1 and b 2 ) and OO 7 (panels c 1 and c 2 ) after 4 (panels a 1 , b 1 and c 1 ) and 12 Myr (panels a 2 , b 2 and c 2 ) since the beginning of subduction. Dashed and dotted lines indicate the gravitational contribution of the sole positive and negative density anomalies, respectively, for the same models in the same configuration. OO 1 : u s =5 cm yr -1 ; θ s =45 • ; c f =1; OO 4 : u s =8 cm yr -1 ; θ s =45 • ; c f =0.5; OO 7 : u s =2 cm yr -1 ; θ s =45 • ; c f =0.25.

Regional gravity pattern
where λ is the ellipsoidal longitude, φ is the ellipsoidal latitude and h is its height above the sea level. Among others, this functional is the most suitable for a direct comparison with our model predictions, because both normal gravity and actual gravity are calculated in the same point, without introducing any other reference point. Fig. 19 includes for each subduction complex three profiles perpendicular to the trench, along which the observed gravity disturbance will be compared to the modelled gravity disturbance. These profiles are chosen in order to sample three different regions of the subduction complex, taken collectively as representatives of the whole arc, being each of them perpendicular to a different sector of the arc itself.
The two settings portray similar gravity features away from the trench region over the upper plate, characterized by a widespread positive gravity disturbance up to 200 mGal in both the continental and oceanic domains. Differences are instead visible over the lower plates, with gravity disturbances very low, close to zero, over the oceanic plate at Sumatra and higher gravity positive values characterizing the upper oceanic plate of the Mariana subduction. In the proximity of the trench, the map shows the bipolar trend typical of the subduction complexes, and beyond the variability depending on the position along each arc, there are some features that are peculiar of each subduction context: ocean-ocean and ocean-continent. For the ocean-continent subduction, the peak-to-peak gravity disturbance varies from approximately 240 mGal in the southern part to 280 mGal in the central portion and to 140 mGal to the north. For the ocean-ocean subduction, the peak-to-peak gravity disturbance is generally higher, approximately 30 per cent, varying from 310 mGal in the south, 270 mGal in the centre and 200 mGal to the north. In addition, for ocean-ocean subduction, the gravity disturbance across the trench is characterized by a well-developed single minimum as low as -220 mGal, and for the ocean-continent subduction, a single minimum is not well defined, and a negative value that is approximately a factor two smaller, at most -120 mGal, is observed. Finally, the lateral extent of the gravitational contribution trough is smaller across the Sumatra (smaller than 200 km) than across the Mariana (larger than 300 km) subduction complex.

Modelled gravity disturbance
In the previous section, we considered the density distribution at the beginning of subduction as the reference density configuration, because the target of the analysis was the perturbation induced by subduction on the gravity field. To properly compare model prediction to the EIGEN-6C4 gravity data, based on a geodetic normal Earth characterized by a uniformly thick continental crust and whose gravity is γ (h, φ) in eq. (8), we here define a 1-D model normal Earth such that it is possible to calculate a modelled gravitational disturbance, which is compliant with the geodetic definition of gravity disturbance given by eq. (8).
Since the normal Earth of the EIGEN6C4 model is based on a continental lithosphere, in order to construct the model normal Earth, we consider the vertical density distribution at a great distance from the trench, on the side of the continent, for all the ocean-continent models. Furthermore, since the normal Earth must represent the present-day reference density distribution, it is built at each time the comparison is made. Fig. 20 shows the density distribution of the model normal Earth obtained after 40 Myr from the beginning of subduction for model OC 3 that will prove to be the best-fitting model for Sumatra.
We thus calculate the modelled gravitational disturbance as the difference between the magnitude of the modelled gravitational contribution of the perturbed masses at a given time t, g model , and the magnitude of the gravitational contribution from the model normal Earth at the same time t, γ model , both calculated at the same point: where x indicates the position of the observational point with respect to the trench and h is its height above sea level. g model is calculated as expressed in eqs (5) and (6), but with ρ e instead of ρ e Among all the gravity functionals, the gravity disturbance is the most suitable one for a direct comparison with our model prediction, because it is based on the difference between the normal gravity and the actual gravity computed in the same point, without introducing any other reference point.

Sumatra subduction
For the Sumatra subduction, we consider all the ocean-continent models with a subduction velocity of 5 cm yr -1 , compatible with the tectonic information of Section 4.1, and we calculate the gravitational contribution of the mass distribution predicted after approximately 40 Myr from the beginning of the subduction, accounting for a 4-km-thick ocean overlying the subducting plate. Below the discussion will be limited to the only model that shows the best agreement with the data, namely, model OC 3 .
Panel a of Fig. 21 compares the gravity disturbance predicted by model OC 3 and the gravity disturbance based on EIGEN-6C4, extracted along three sections perpendicularly crossing the Sumatra trench in the north (dashed red line), in the centre (solid red line) and in the south (dotted red line). The model (solid black curve) based on a coupling factor of c f =0.5 and on the bathymetry of the southern profile (Fig. 21b) fits the envelope of the profiles well in terms of positioning and height of the maxima of 30 and 200 mGal (solid red curve), with a depth of the trough or gravitational disturbance minimum of -120 mGal (dashed red curve) and with the asymptotic average values of a few tens of mGal over the subducting and overriding plates, negative and positive, respectively. As expected for a model that does not account for processes except subduction, the series of local maxima are not reproduced. Instead, the model does reproduce the narrowest trough of approximately 175-180 km of the red solid curve, evaluated from the distance between the two positive peaks surrounding the trough, to be compared in the following with that of the ocean-ocean subduction.
Panel (c) of Fig. 21 depicting the marker distribution and the density anomalies after 40 Myr since the beginning of subduction, provides the explanation for this gravitational disturbance pattern due to the narrow region of negative density contrast, as large as 80 km at most, to a depth of approximately 80 km, embedded between the subduction cold plate and the upwelled lithospheric mantle. This mantle substitutes the less dense continental crust, as indicated by the white region below the thinned crust, dark and light red in panel c of the marker distribution between 20 and 40 km, thus producing a positive density contrast, portrayed by the shallow red zone in the inset of panel (c). This shallow positive density contrast is thus effective in the thinning of the trough in the ocean-continent subduction. Furthermore, this positive density contrast has the effect of increasing the positive gravitational peak located on the continental side. This positive density contrast at depths of 20-50 km originating from the indentation of dense upper mantle material within the continental crust in this ocean-continent case to some extent shields the deeper negative density contrasts, whose contribution to the formation of the trough is thus inhibited.

Mariana subduction
For the Mariana subduction, we consider all the ocean-ocean models with the same value of subduction velocity of 5 cm yr -1 as for Sumatra, and we calculate the gravitational contribution of the mass distribution after 45 Myr, compatible with Sumatra, accounting for 5-km-thick ocean overlying the subducting plate and a 4-km-thick ocean overlying the overriding plate. Below the discussion will be limited to the only model that shows the best agreement with the data, namely, model OO 1 .  Fig. 22 shows the same information as Fig. 21 but for the Mariana subduction. A wider trough with respect to the ocean-continent subduction characterizes this region and is well reproduced by the OO 1 model. This model accounts for a coupling factor of c f =1 and builds on the bathymetry of the southern profile (Fig. 22, panel b). With respect to the ocean-continent subduction, the ocean-ocean subduction is also characterized by a large recirculation of light material in the mantle wedge, which is responsible for the wide region of negative density contrasts in blue in the inset of panel (c) of Fig. 22, encompassing a region in the horizontal direction as large as approximately 150 km and leading to the prediction of a large gravity trough in agreement with the data (Fig. 22a). Another major difference that makes the distribution of light material more effective in the ocean-ocean subduction in the formation of a pronounced and wide gravitational trough relies on the fact that in this type of subduction, in contrast to the ocean-continent one, the dense mantle does not substitute the lighter continental crust, which causes the formation of the positive density contrast observed above, contributing to the narrowness of the trough in the OC 3 model shown in Fig. 21.
Panel ( Figure, we can provide an explanation for other characteristics of the ocean-ocean case, namely, the second-order maxima characterizing the gravitational signature within the trough on the side of the overriding plate, both along two data profiles and along the model black curve. The recirculation of light crustal material and sediments indicated by the markers in panel (c) of Fig. 22 is characterized by a core of small-scale positive density contrast originating from the dense oceanic subducting crust, carried at shallow depth by the vigorous convective circulation occurring in the mantle wedge once out-scraped from the top of the slab. Our results allow us to provide a physical explanation for the occurrence of the small-scale local maxima in the gravitational disturbance on the side of the subducting slab due to localized dense material from the out-scraping of the slab.

C O N C L U S I O N S
We have developed a thermomechanical model that contains most of the complexities, such as compositional and phase changes, hydration, sedimentation and the degree of plate coupling, within the frame of a convective viscoplastic model, to analyse the static and time -dependent gravitational contribution of subduction. Each of these complexities plays a role in shaping the pattern of the gravitational disturbance, particularly in the trench region, when a self-consistent comparison with EIGEN-6C4 model is also performed. The effects on the gravity pattern due to a variation of several parameters (such as the degree of plate coupling, the subduction velocity, the subduction dip and the ocean-ocean or ocean-continent environment) have been also exploited. The model results indicate that: (i) As the plate coupling decreases, (a) the cooling of the wedge area intensifies while, at great depths, the subducted slab remains warmer; (b) the horizontal extension of the hydrated wedge area increases as the plate coupling decreases, while the maximum depth decreases; (c) the shallow negative density anomaly, associated with the trench trough, the subducted and recycled crust and sediments and the mantle hydration, increases and widens; (d) the gravitational contribution from the negative density anomalies increases as well as the amplitude of the gravity trough; (ii) As the subduction dip angle increases: (a) the p-T condition becomes unfavorable for mantle hydration and recirculation of light material; (b) intense but localized shallow negative and positive density anomalies develop; (c) more intense shallow positive density anomalies while less intense deep positive density anomalies develop; (d) the gravitational contribution from the negative density anomalies decreases as well as the amplitude of the gravity trough; (e) the gravitational contribution from the positive density anomalies increases; (iii) As the subduction velocity decreases (a) a thermal thickening of the subducted slab is enhanced at shallow depths; (b) a widespread positive density anomaly develops at the base of the upper plate and along the subducted slab; (c) the negative density anomaly is not affected by changes in the subduction velocity; (d) the gravitational contribution due to the positive density anomalies increases as well as the width of the gravity trough. When the rate of changes are normalized with respect to the corresponding subduction velocity, non-significant differences are detectable and it is possible to estimate a common rate of change of 0.008 μGal yr -1 (cm yr -1 ) -1 ; (iv) The differences between ocean-ocean and ocean-continent environments are due to the different dynamics within the wedge embedded between the subducting and overriding plates for the two types of subduction, affecting in particular the amount of low-density material and thus the gravity low characterizing the trough. The gravity trough is generated by the negative topography and by the low-density material that overcomes, in proximity to the trench, the global positive gravitational disturbance caused by the cold subducted lithosphere, which remains unaffected by the nature of the overriding plate, being a continent or an ocean; (v) A gravity change of approximately 60 mGal occurs within about 2 Myr when the subducted slab crosses the olivine-spinel transition.
The study of the EIGEN-6C4 gravitational disturbance patterns of the Sumatra and Mariana subductions have allowed us to strengthen the analysis of the gravitational signature in ocean-continent and ocean-ocean subductions in terms of the physics of the processes occurring during the convergence of the plates. Although both types of subduction show the classical 2-D dipolar profile perpendicular to the trench, they can be distinguished in terms of some fundamental features.
Our modelling is able to reproduce the gravity disturbance difference of 250-300 mGal well between the maximum and the minimum, characterizing both types of subduction. In the same way it reproduces the fundamental differences highlighted by the EIGEN-6C4 data: (a) the width of the trough (larger for the ocean-ocean subduction than for the ocean-continent one); (b) the symmetry, in terms of the different amplitudes of the two positive gravity peaks facing the trench.
Although focused on the analysis of the different gravitational disturbance patterns of the two types of subduction, our study provides a physical explanation for the broadness of the negative gravitational contribution for mature subductions (as the Mariana) compared to immature ones, as highlighted by Kim et al. (2009). Serpentinization, partial melt due to frictional heating in the mantle wedge and differences in the density structure between the overriding and subducting plates were in fact proposed by Kim et al. (2009) to explain the broad gravitational contribution in the trench region for mature ocean-ocean subductions.
Finally, our results show that gravity patterns and their rate of change in subduction zones provide important information not only on their anomalous density structure but also on the dynamics of the subduction process, complementing those studies in which density anomalies are obtained from seismic tomography or petrological information (e.g. Kim et al. 2009;Simmons et al. 2010;Moulik & Ekström 2016).

A C K N O W L E D G E M E N T S
The work was supported by the ESA funded project Gravitational Seismology ITT AO/1-9101/17/INB. The authors thank the Editor J. C. Afonso and the three reviewers, A. Negredo and two anonymous. All figures have been made using GMT -The Generic Mapping Tools (Wessel et al. 2013).

A P P E N D I X A : M AT H E M AT I C A L F O R M U L AT I O N
The eqs (1), (2) and (3) are numerically integrated via the 2-D FE thermomechanical code SubMar (Marotta et al. 2006), which uses the Petrov-Galerkin method to integrate the conservation of energy equation and the penalty formulation: to solve for pressure, where λ is the penalty parameter. We follow Sewell (1981)'s approach and fix λ = μ √ r p , where r p denotes the machine relative precision and μ is interpreted as a bulk viscosity (Hughes et al. 1979;Donea & Huerta 2003;Marotta et al. 2006;Thieulot 2014). The marker-in-cell technique (e.g., Christensen 1992) has been used to compositionally differentiate the different types of materials. At the beginning of the evolution, 455 489 markers identified by different code are spatially distributed, with a density of 1 marker per 0.25 km 2 , to define the upper and lower oceanic crust, the upper and lower continental crust and the continental mantle. During the evolution of the system, each particle (marker) is advected using a fourth-order (both in time and in space) Runge-Kutta scheme. At each time, the elemental density of each type of particle defines the composition of each element of the grid, C e , such as where where N e i and N e 0 are the number of particles of type i within element e and the maximum number of particles of any type that the element e contains, respectively. Each elemental property P e (conductivity, specific heat at a constant pressure, density and viscosity) depends on the composition such that within each element e, P e can be expressed as follows: where P m and P i are the property for mantle and for any particles of type i, respectively. For what concerns the density, if it depends only on temperature, the corresponding variations are generally small, and the fluid can be assumed to be incompressible. The density can consequently be treated as a constant in the continuity (1) and energy eq. (3), leading to an incompressible fluid (e.g. Christensen & Yuen 1985;Ismail-Zadeh & Tackley 2010). If density variations depend not only on temperature variations but also on pressure variations as considered herein, then compressibility must be taken into account. Accounting for compressibility can be achieved by using either the extended Boussinesq approximation or the anelastic approximation (Gerya 2010;Ismail-Zadeh & Tackley 2010). Here, we implemented the extended Boussinesq approximation in the following form: where ρ m (T, p) is the mantle density; ρ i (T, p), with i varying from 1 to the total number of no-mantle material types included in the model, are the densities of no-mantle material (e.g., Ismail-Zadeh & Tackley 2010;Thielmann & Kaus 2012;Quinquis & Buiter 2014;Thieulot 2014); and C i are the corresponding compositions. Both ρ m and ρ i are functions of temperature, pressure and material type. The extended Boussinesq approximation allows us to consider the density as a constant in the continuity eq. (1), while it accounts for density variations in the buoyancy term of the momentum eq. (2) and for both the shear and adiabatic heating terms in the energy eq. (3), where the Table A1. Values of the material and rheological parameters used in the analysis. References: (a) Ranalli & Murphy (1987); (b) Afonso & Ranalli (2004); (c) Kirby (1983); (d) Haenel et al. (1988); (e) Chopra & Peterson (1981); (f) Dubois & Diament (1997) and Best & Christiansen (2001); (g) Roda et al. (2011); (h) Schmidt & Poli (1998); (i) Gerya & Stockhert (2006); (j) Roda et al. (2012); (k) Gerya & Yuen (2003); (l) Borghesan (2011); (m) Schmeling et al. (2008 where H r is the radiogenic heating, H s = τ i j˙ is the shear heating and H a = T α D p Dt , where α is the thermal expansion coefficient, is the adiabatic heating per mass unit (e.g. Christensen & Yuen 1985;Gerya 2010;Ismail-Zadeh & Tackley 2010;Thielmann & Kaus 2012).
The extended Boussinesq approximation is widely used when dealing with compressible fluids (Gerya et al. 2004;Burg & Gerya 2005;Gorczyk et al. 2006;Duretz et al. 2011;Wang et al. 2019) accounting for changes in the buoyancy term due to density changes in the momentum and in the energy equations, including those from phase transformations, and for a time-invariant density in the continuity equation (Ismail-Zadeh & Tackley 2010). In order to enforce locally and only for lithospheric material the mass conservation where phase changes occur, Afonso & Zlotnik (2011) introduced in the momentum equation the contribution from the divergence of the velocity field obtained by approximating the continuity equation according to ∂u/∂x + ∂v/∂y −(1/ρ)( ρ/ t).
ρ denotes in each Lagrangian particle the difference between the density in the current and previous time step, which generates a localized velocity field that compresses or expands the material at the time of the phase transition in the lithosphere, enforcing locally the mass conservation and providing a higher order approximation within the scheme of the extended Bousinnesq approximation. We also account for the changes in time of the density due to phase changes, by modifying at each time step the current density according to the pressure and temperature conditions, as required by the extended Bousinessq approximation, both in the buoyancy term of the momentum equation and in the energy equations, as described in Appendix C. In our model mass conservation is guaranteed in the study domain by the horizontal material flow in and out of the domain across its vertical left boundary, as shown in Fig. 2. The model combines a linear viscous rheology for the sublithospheric mantle with a linear viscoplastic rheology for the lithosphere. To this purpose, for each material type i, we compute the following two viscosities: where μ 0, i and E i are the reference viscosity at the reference temperature T 0 and the activation energy, respectively, for material type i with μ 1 , μ 2 and μ 3 computed combining the Byerlee's law criterion and the Tresca criterion as implemented in Regorda et al. (2017) where p is the pressure,˙ 1 ,˙ 2 and˙ 3 are the principal strain rates and σ Y is the yield stress, defined by the simplified formulation of Byerlee's law criterion, σ Y = β · y, with β = 40 MPa km −1 and y depth in kilometres. Outside the lithosphere, the effective viscosity is then calculated as follows: Within the lithospheric layer the effective viscosity is instead defined as follows: Table A1 lists the values of the material and rheological parameters used in the analysis. Instantaneous erosion/sedimentation mechanism is simulated as in Roda et al. (2010) by using the substitution technique, in which all the crustal particles lying above a prescribed height (h t in Table 1) are replaced with air particles and, at the same time, an equal number of water particles buried into the trench region are transformed into sediments. This procedure reproduces an erosion and sedimentation rate variable in time, as function of the system dynamics. Further details can be found in

A P P E N D I X B : H Y D R AT I O N A N D S E R P E N T I N I Z AT I O N O F T H E W E D G E A R E A
Water contained in hydrous phases in H 2 O-saturated MORB basalt can be transported to large depths (e.g. Liu et al. 2007;Faccenda 2014;Rosas et al. 2016), where the oceanic plate is then dehydrated because of the increased temperature and pressure (Schmidt & Poli 1998;Liu et al. 2007;Faccenda et al. 2009;Faccenda & Mancktelow 2010). Schmidt & Poli (1998) and Liu et al. (2007) suggested that complete dehydration is achieved at a depth between 70 and 350 km. Furthermore, Faccenda et al. (2009), Faccenda & Mancktelow (2010 and Faccenda (2014) showed that tectonic stresses can influence the hydration pattern in the subducted slab, allowing fluids to penetrate into the slab and favouring their transport to great depths up to the base of the upper mantle. The presence of hydrothermal circulation in the oceanic crust can affect the thermal state in the subducting slab, which generates higher hydrous phase stability at higher depths compared with those predicted by classical thermal models of the subduction zone (Perry et al. 2016;Rosas et al. 2016). However, hydrothermal circulation does not have a significant impact on the depth of dehydration reactions at shallow depths (Rosas et al. 2016).
The phases with large contributions to the water budget of subducting MORB are lawsonite, chlorite and amphibole; however, at pressures higher than 3-4 GPa, lawsonite constitutes the only hydrous phase. At pressures above the zoisite stability field, the lawsonite breakdown reaction has a positive dp/dT slope in the coesite stability field, while it has a negative dp/dT slope in the stishovite stability field, above 8 GPa (Schmidt & Poli 1998). In our models, the maximum dehydration depth of the oceanic crust, below which the water content in the subducting oceanic plate is negligible, has been determined using the stability field of lawsonite as follows: where p dehydr is the maximum pressure at which lawsonite is stable calculated for each oceanic crustal marker with temperature T marc . The water released by the slab hydrates the mantle wedge and can lead to its serpentinization, with a consequent decrease in viscosity and density (Honda & Saito 2003;Gerya et al. 2002;Arcay et al. 2005). The rheological weakening of the mantle wedge has been simulated by assuming a constant viscosity of 10 19 Pa · s and a density of 3000 kg m -3 for the serpentinized mantle (Honda & Saito 2003;Arcay et al. 2005;Gerya & Stockhert 2006;Roda et al. 2010).
In our models, the stability field of serpentine is calculated for each element as follows: where p hydr represents the maximum pressure at which serpentine is stable at the elemental temperature T elem . The hydrated area is limited from below by the oceanic subducting plate. Specifically, to better delineate the geometry of the lower border of the hydrated area, the subducted plate is subdivided into segments of equal length, and the deepest dehydrated oceanic crust marker for each segment is identified at each time during the system's dynamic evolution. The line that connects these markers defines the lower limit of the hydrated area.
We implement a free-water migration mechanism at a constant velocity of 20 cm yr -1 , in accordance with the velocities obtained by Quinquis & Buiter (2014) for different water migration schemes (10-70 cm yr -1 ). Further, Quinquis & Buiter (2014) demonstrated that the exact manner of water migration does not impact the dynamics in the mantle wedge. We also assume that the used time step (50 000 yr) is sufficient to have a complete serpentinization of the hydrated areas (Arcay et al. 2005).

A P P E N D I X C : P H A S E C H A N G E S
In all models, the phase changes were introduced both for the mantle and for all the lithologies considered for the crust. For each phase change, we introduced density variations in relation to the variation in pressure and temperature (p-T) conditions predicted by the models.
Regarding the phase changes occurring in the mantle, Fig. C1 shows the assumed stability fields for olivine, wadsleyite, ringwoodite and perovskite+periclase. Considering a temperature of ∼1600 K for the upper mantle, olivine is expected to transform into wadsleyite at a pressure of ∼13 GPa, while the wadsleyite → ringwoodite transformation is expected to occur at ∼17.5 GPa (Katsura & Ito 1989;Kerschhofer et al. 1998), both with a positive Clayperon slope. We considered a density for the olivine of 3200 kg m -3 (Deer et al. 1992), while the α → β and β → γ phase transformations result in increases in density of 6 and 2 per cent (Kerschhofer et al. 1998), respectively, obtaining densities of 3392 and 3460 kg m -3 for the wadsleyite and ringwoodite, respectively. At the bottom of the upper mantle, ringwoodite transforms into perovskite + periclase (MgSiO 3 +MgO). At a temperature of 1600 K, this transformation occurs at ∼22 GPa (corresponding to the 660-km seismic velocity discontinuity in the mantle), with a negative Clayperon slope (Fei et al. 2004) and refs. therein), and produces an additional increase in density. In particular, considering densities of ∼4000 and ∼3600 kg m -3 for the perovskite and the periclase (Deer et al. 1992), respectively, we assumed a density of 3800 kg m -3 for the lower mantle.
In contrast, Perple X (Connolly 1990) has been used to introduce density variations related to multiple phase changes occurring at the same time in the oceanic crust and continental crust and in the sediments, each characterized by particular bulk compositions (Borghesan 2011): (i) the oceanic crust has the bulk composition of a mid-Atlantic ridge N-MORB (Schilling et al. 1983); (ii) the continental crust has been separated into the upper crust, with the bulk composition of the Mucrone granite (Oberhansli et al. 1985), and the lower crust, with the bulk composition of the Ivrea diorite (Bigioggero et al. 1979); (iii) the sediments have the bulk composition of the Tre Valli Bresciane micaschist (Origoni & Gregnanin 1983).
For each bulk composition, Perple X produced 10 000 different densities for temperatures between 300 and 1600 K and pressures between 0 and 7 GPa (Borghesan 2011 and Fig. C2), and each density has been linked to different p-T conditions predicted by the models. In this way, the variations in density inside the oceanic crust, the continental crust and the sediments are nearly continuous in relation to variations in the p-T conditions (Fig. C2). At each time step we iterate between the density and the pressure, driving the density changes due to phase change.
Phase changes have a dual effect, determining a variation in density and release or absorption of latent heat. However, in our models, we did not include in eq. (3) the latent heat produced by phase changes occurring in the continental and oceanic crust, by phase changes related to hydration of the mantle wedge and by phase changes regarding the transformation of the olivine in the deep mantle. Afonso & Zlotnik (2011) demonstrated that the impact of the crustal metamorphic reactions on the local thermal state is on the order of 5-10 • C, and then, the latent heat produced by phase changes of the crust can be ignored in large-scale models. Regarding the phase changes in the mantle wedge, the latent heat of serpentinization can be neglected because it is balanced at large scales by the latent heat produced by de-serpentinization (Perez-Gussinye & Reston 2001; Rupke et al. 2004Rupke et al. , 2013. In addition, the influence of latent heating should laterally impact the dehydration fronts, with a slight impact on the width and on the mechanisms related to the hydrated area, such as the erosion of the upper plate and the activation of a small-scale convective flow (Arcay et al. 2005;Yamato et al. 2007). The main phase transformations occurring in the upper mantle regard the olivine (Mg, Fe) 2 SiO 4 , which transforms from the α − phase to the high-pressure polymorphs wadsleyite (β − phase) and ringwoodite (γ − phase). These reactions, such as the olivine-wadsleyite phase transition, can increase the local temperature up to 60-70 • C ( van Hunen et al. 2001;Afonso & Zlotnik 2011); however, Negredo et al. (2004) observed that the introduction of the latent heat in thermal models determines a variation in the maximum depth reach by the α − phase of less than 50 km inside the slab, while no differences have been observed outside the slab. Moreover, Blom (2016) observed that the latent heat produced at the olivine → wadsleyite transformation is one order of magnitude smaller than the shear heating and two orders of magnitude smaller than the adiabatic heating. Thus, the effects of the introduction of the latent heat can be considered negligible for what concerns the large-scale dynamics of the upper mantle and the Bouguer gravity anomalies corresponding to the subduction zone. Figure C1. Stability fields of olivine (α), wadsleyite (β), ringwoodite (γ ) and perovskite+periclase in isochemical peridotitic mantle. Based on Katsura & Ito (1989) and Fei et al. (2004).

A P P E N D I X D : C O U P L I N G FA C T O R T E C H N I Q U E
To implement the degree of coupling, we started from the classic split node technique, which allows for each node belonging to the discontinuity plane to be assigned of two different velocities, depending on the element to which it belongs. The split node technique, originally introduced by Jungels (1973) and Jungels & Frazier (1973) for elastic rheology, is a method used to introduce a discontinuity in the displacement/velocity along fault planes within numerical computations based on the FE method. Normally, the displacement/velocity at nodes shared by more than one element is the same for each element. The split nodes are special nodes in which the displacement/velocity is different depending on the element it belongs to. In this way, a generic differential slip/slip-rate + u can be assigned to a node when referring to the element on one side of the fault, while a differential slip/slip-rate − u is assigned to the same node when referring to the element on the other side of the faults (yellow circles in panel b of Fig. 2). The fault is thus defined by a line of split nodes. This technique is implemented locally, thus simplifying the introduction of split nodes and the assembly of the vector of the charges and the stiffness matrix, and no degree of freedom is added to the nodes (Melosh & Raefsky 1981). The original formulation of the split node technique, however, has not been found suitable for our modelling because it requires that the slip/slip-rate is known a priori. In our case, it is not possible to know a priori the amount of slip rate along the fault plane, leading to the results from the self-consistent thermomechanical evolution of the model.
We have modified the original formulation of the split node technique by introducing a coupling factor c f that varies from 0 to 1 and indicates the percentage difference between the velocities at the two split nodes in such a way that where u l and u r are the velocities of the same split node belonging to the left and right elements, respectively (Fig. 2). Unlike the classical split node technique, which only involves modifying the load vector, this new technique also requires modification of the stiffness matrix.
Considering the two quadratic triangular elements of Fig. D1 and assuming that 1 of the 6 nodes is a split node, namely, node k (red circle in Fig. D1) with a coupling factor c f (k) different from 1. The velocity calculated in node k is equal to u(k) when node k is considered to belong to the left element e l , while it is equal to c f (k) · u(k) when node k is considered to belong to the right element e r . To simplify the implementation of the new technique, we define a coupling factor to all the nodes of the numerical grid, and we initialize the array with the default value equal to 1. For the slip node, the value of the coupling factor changes into the actual value only when the slip node is treated as belonging to the right element. Its value remains unchanged when the same slip node is treated as belonging to the left element. Following the standard formulation of the FE method, at the local scale of each element, the expression K e · u e = F e becomes: ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ K 11 K 12 K 13 K 14 K 15 K 16 K 21 K 22 K 23 K 24 K 25 K 26 K 31 K 32 K 33 K 34 K 35 K 36 K 41 K 42 K 43 K 44 K 45 K 46 K 51 K 52 K 53 K 54 K 55 K 56 where K e is the element stiffness matrix, u e is the vector of the element nodal displacement and F e is the vector of the element generalized nodal forces.