Diane Arcay
email: diane.arcay@gm.univ-montp2.fr

Modelling the interplate domain in thermo-mechanical simulations of subduction: Critical effects of resolution and rheology, and consequences on wet mantle melting

Keywords: subduction interface, thermo-mechanical modelling, rheology, numerical resolution, oceanic crust buoyancy, wet mantle partial melting

The present study aims at better deciphering the different mechanisms involved in the functioning of the subduction interplate. A 2D thermo-mechanical model is used to simulate a subduction channel, made of oceanic crust, free to evolve. Convergence at constant rate is imposed under a 100 km thick upper plate. Pseudo-brittle and non-Newtonian behaviours are modelled. The influence of the strength contrast modelling the subduction channel, represented by the difference in activation energy between crust and mantle (∆E a ) is investigated to examine in detail the variations in depth of the subduction plane down-dip extent, z coup . First, simulations show that numerical resolution may be responsible for an artificial and significant shallowing of z coup if the weak crustal layer is not correctly discretized. Second, if the age of the subducting plate is 100 Myr, subduction occurs for any strength reduction ∆E a . The stiffer the crust is, the shallower z coup is (60 km depth if ∆E a = 20 kJ/mol) and the hotter the fore-arc base is. Conversely, imposing a very weak subduction channel (∆E a > 135 kJ/mol) leads to an extreme mantle wedge cooling and inhibits mantle melting in wet conditions. Partial kinematic coupling at the fore-arc base occurs if ∆E a = 145 kJ/mol. If the incoming plate is 20 Myr old, subduction can occur if the crust is either stiff and denser than the mantle, or weak and buoyant. In the latter case cold crust plumes rise from the subduction channel and ascend through the upper lithosphere, triggering (1) partial kinematic coupling at the fore-arc base, (2) fore-arc lithosphere cooling, and (3) partial or complete hindrance of wet mantle melting. z coup then ranges from 50 to more than 250 km depth and is time-dependent if crust plumes form. Finally, subduction plane dynamics is intimately linked to the regime of subduction-induced corner flow. Two different intervals of ∆E a are underlined: 80-120 kJ/mol to reproduce the range of slab surface temperature inferred from geothermometry, and 10-40 kJ/mol to reproduce the shallow hot mantle wedge core inferred

Introduction

The mantle wedge tip and subduction plane are key structures of the subduction factory, affecting not only the subduction viability but also seismogenic properties of the subduction plane, fluid transfers, and melt migration towards the volcanic arc (Fig. 1). Using a thermo-kinematic model applied in 17 subduction zones constrained by heat flux profiles, [START_REF] Wada | Common depth of slab-mantle decoupling: Reconciling diversity and uniformity of subduction zones[END_REF] suggested that the Maximum Depth of kinematic Decoupling (MDD) along the subduction plane was constant for all considered subduction zones and equal to 80 km depth, the process controling it remaining to be determined. Other observations advocate on the contrary for a certain degree of variability, as indicated for instance by the depth to slab under volcanic front ("sub-arc slab depth", Fig. 1) : 110±20km in [START_REF] Jarrard | Relations among subduction parameters[END_REF], varying more broadly in recent studies: 60-120 km depth in [START_REF] England | A simple analytical approximation to the temperature structure in subduction zones[END_REF] and 70-170 km in [START_REF] Syracuse | Global compilation of variations in slab depth beneath arc volcanoes and implications[END_REF]. Even if the average sub-arc slab depth is always centered at 105 km, corresponding to the most frequent overlapping plate thickness [START_REF] Heuret | Plate motions, slab dynamics and back-arc deformation[END_REF], the volcanic arc location may vary as a function of convergence parameters, even if the exact governing process is debated [START_REF] Grove | Kinematic variables and water transport control the formation and location of arc volcanoes[END_REF][START_REF] England | Melting above the anhydrous solidus controls the location of volcanic arcs[END_REF]. To assess the exact MDD, heat flow data are useful: from a minimum at trench, heat flux increases landwards to a maximum plateau assumed to start right above the MDD [START_REF] Honda | Thermal structure beneath Tohoku, northeast Japan-a case study for understanding the detailed thermal structure of the ubduction zone[END_REF][START_REF] Furukawa | Depth of the decoupling plate interface and thermal structure under arcs[END_REF]. Nonetheless, the assertion by [START_REF] Wada | Common depth of slab-mantle decoupling: Reconciling diversity and uniformity of subduction zones[END_REF] of a MDD uniformly equal to 80 km depth is based on heat-flow data presenting either an extremely low spatial density or a significant scattering. Even if the down-dip extent of subduction planes might actually be located at 80 km depth for all subduction zones, the observations used to state it are not completely convincing yet, whereas the variability of well-constrained sub-arc slab depths is meaningful. In any case, the processes regulating the subduction plane geometry have to be better understood. Notable progress have been made in understanding the subduction plane dynamics. On one hand, van [START_REF] Van Keken | High-resolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle[END_REF] using thermo-kinematic models showed that modelling a temperature-and stress-dependent viscosity drastically increased thermal gradients at the slab surface. [START_REF] Wada | Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization[END_REF] and [START_REF] Hirauchi | Rheological contrast between serpentine species and implications for slab-mantle wedge decoupling[END_REF] showed that the formation of low-temperature hydrated minerals such as serpentine in the thermal boundary layer covering the subducting slab reduces viscous couplings between slab and asthenosphere and cools the mantle wedge. Thermo-kinematic models have been pointed out by geochemists studying mantle melting at the slab surface (e.g., [START_REF] Kelemen | Thermal structure due to solid-state flow in the mantle wedge beneath arcs[END_REF][START_REF] Plank | Emerging geothermometers for estimating slab surface temperatures[END_REF][START_REF] Grove | The role of H 2 O in subduction zone magmatism[END_REF]Cooper et al., 2012, Fig. 1), this because simulated slab surfaces are hotter than obtained in thermo-mechanical models, in better agreement with geochemical analyses. The cause likely comes from the prescribed rigid lid within the upper plate that partly includes the underneath low viscosity wedge, which should otherwise flow onto the slab surface (Arcay et al., 2007b;[START_REF] Wada | Dynamics of subducting slabs: Numerical modeling and constraints from seismology, geoid, topography, geochemistry, and petrology[END_REF], and cool it. [START_REF] Conder | A case for hot slab surface temperatures in numerical viscous flow models of subduction zones with an improved fault zone parameterization[END_REF] proposed a parameterization of the subduction fault taking into account temperature-and strain-rate dependence to limit the arbitrariness of the imposed rigid lid thickness. Moreover, partial kinematic coupling has been assumed to occur at the fore-arc base [START_REF] Kneller | B-type olivine fabric in the mantle wedge: Insights from high-resolution non-Newtonian subduction zone models[END_REF][START_REF] Kneller | Stress, strain, and B-type olivine fabric in the fore-arc mantle: Sensitivity tests using high-resolution steady-state subduction zone models[END_REF][START_REF] Van Keken | Subduction factory: 4. depth-dependent flux of h2o from subducting slabs worldwide[END_REF], i.e., velocities in the thin layer covering the down-dip extent of the subduction channel would be much lower than corner flow velocity, typically about a few percents of subduction speed (Fig. 1). The origin of partial kinematic coupling being unclear, a possible explanation is proposed in this paper. On the other hand, thermo-mechanical models showed that if the subduction plane was free to evolve its structure (depth of brittle-ductile transition and MDD) depended on the modelled brittle rheology and on ductile strength, as well as on subduction speed and age at trench of the subducting plate [START_REF] Arcay | Overriding plate thinning in subduction zones: Localized convection induced by slab dehydration[END_REF][START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF]. However, the influence of interplate rheology, addressed by [START_REF] Wada | Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization[END_REF] by imposing the MDD, has not been completely established in a self-consistent model, although interactions between interplate shearing and corner flow dynamics may be crucial to regulate the subduction interplate. This motivates the present study. In thermo-mechanical models, the tangential velocity discontinuity along the subduction interface (Fig. 1) is classically modelled by a thin layer weaker than the surrounding mantle by several orders of magnitude (e.g., [START_REF] Kopitzke | Finite element convection models: comparison of shallow and deep mantle convection, and temperatures in the mantle[END_REF][START_REF] Gurnis | Controls on the structure of subducted slabs[END_REF][START_REF] King | The relationship between plate velocity and trench viscosity in Newtonian and power-law subduction calculations[END_REF][START_REF] Gerya | Why is terrestrial subduction one-sided?[END_REF][START_REF] Quinquis | The role of boundary conditions in numerical models of subduction dynamics[END_REF][START_REF] Agrusta | The effect of metastable pyroxene on the slab dynamcis[END_REF][START_REF] Garel | Interaction of subducted slabs with the mantle transition-zone: A regime diagram from 2-d thermo-mechanical models with a mobile trench and an overriding plate[END_REF][START_REF] Crameri | Parameters controlling dynamically selfconsistent plate tectonics and single-sided subduction in global models of mantle convection[END_REF][START_REF] Marques | Speculations on the impact of catastrophic subduction initiation on the earth system[END_REF]. A similar set-up is used in this paper to study carefully how the subduction channel geometry is regulated and how its strength affects the global subduction zone structure. As the subduction plane is likely to equilibrate at the time-scale of tens of Myr, long-term experiments are performed (∼20 to 40 Myr of convergence). If the upper plate is thick (∼100 km), forming a viable subduction channel is tough when the incoming plate is young and buoyant [START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF]. That is why the oceanic crust density will also have to be tested when modelling a young plate subduction. Results show that modifying both the subduction channel strength and density generates a broad variety of subduction regimes. At last, as numerical resolution can be critical [START_REF] Schmeling | A benchmark comparison of spontaneous subduction models -Tonwards a free surface[END_REF], it has to be investigated. These different reasons finally yield a long study but, I hope, consistent. The paper starts with a short summary of key processes implied in the functioning of a subduction interplate. Then, the model set-up is described. Next, the strong influence of numerical resolution is detailed first. Afterwards, the paper presents what kind of thermal structure is obtained in the vicinity of the subduction plane down-dip extent as a function of rheology for two endmember subducting plate ages at trench (either 100 Myr, section 5.1, or 20 Myr, section 5.2). The case of an middle age incoming plate (60 Myr) is briefly described to derive a comprehensive regime diagram at the end of the result section. Three types of outputs are described: (1) depth of maximum decoupling along the interplate and velocity field in the mantle wedge, (2) subduction viability, and (3) prediction of wet mantle melting. This study shows that the subduction plane structure and extent are intimately linked to the regime of the subduction-induced corner flow in the asthenosphere. In the last part of the paper, by combining modelling outputs with two main observations typically used to discuss the mantle wedge thermal state (pressure-temperature conditions at the slab surface and spatial extent of the cold fore-arc nose), the range of realistic ductile strength to impose along the subduction plane will be discussed. To conclude, the study argues for an interplate decoupling depth dependent on the subduction setting and likely to vary from a subduction zone to another.

Subduction interplate functioning inferred from previous studies

To sustain subduction, the interplate channel has to allow simultaneously for two opposite processes [START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF]. On the one hand, the velocity discontinuity at shallow depth is maintained by a highly localized shearing (Fig. 1). Strain localization depends on the strength contrast between the soft crustal layer and the surrounding mantle. On the other hand, next to the upper plate base, localized strain has to jump from the subduction channel to the asthenosphere to thermally decouple converging plates [START_REF] Doin | Subduction initiation and continental crust recycling: the roles of rheology and eclogitization[END_REF]. This occurs if the subduction channel is stiffer than hot asthenosphere. Reversely, if the subduction channel is not significantly weak, the asthenospheric flow focussing at mantle wedge tip is needed to mechanically separate converging plates. This phenomenon was well exemplified by [START_REF] Eberle | A numerical study of the interaction between the mantle wedge, the subducting slab, and overriding plate[END_REF]. If the subduction fault is modeled solely by a kinematic boundary condition, and not mechanically (such as by low stress or viscosity condition), when a no-slip condition is imposed on top surface, then the asthenosphere rises along the subduction plane up to z ∼ 0. Indeed, the only means for the flow to satisfy the two kinematic conditions is to propel the weak asthenosphere along the kinematic fault, up to the surface. In this case the asthenosphere acts as the missing decoupling agent. On the contrary, a very weak subduction channel, and/or the absence of corner flow formation prevent localized strain from jumping towards the asthenospheric wedge tip. As a result, the continuous localization of convergence within the weak interplate channel yields horizontal subduction, as modelled for a young and rather buoyant subducting lithosphere [START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF].

Model set-up

Simulations are performed in 2D by solving the equations of mass, momentum, and energy conservation, using the thermo-chemico-mechanical code of [START_REF] Christensen | An Eulerian technique for thermomechanical modeling[END_REF], based on spline-finite elements on an Eulerian grid. Fluid incompressibility is assumed except in the momentum equation. A uniform volumetric heating, as well as shear heating and adiabatic heating are included in the energy equation (extended Boussinesq approximation). Compositional tracers, with specific density and rheology, are advected as a function of the velocity field (fourth order Runge-Kutta scheme, van [START_REF] Van Keken | A comparison of methods for the modeling of thermochemical convection[END_REF].

Boundary conditions and initial state

The 2D simulation box, 2220 km wide and 555 km deep, is filled with two different lithologies: a ∼7 km thick layer overlies the incoming lithosphere, while the rest of the box is made of mantle (Fig. 2). The initial thermal structure of the overriding lithosphere results from a preliminary experiment in which thermal equilibrium between heating from below and cooling from above was performed without plate kinematics, but keeping identical other thermo-mechanical parameters. For simplicity, the upper lithosphere structure is the same for all simulations. A constant convergence rate of 6.5 cm/yr is applied on the top surface far away from the subduction zone (Fig. 2). The crust layer covering the incoming plate surface also fills the initial subduction plane whose dip angle is 30 o until 55 km depth. As crustal rocks are assumed to be significantly weaker than mantle (see below), strain will focus within the initial crustal interplate at convergence start (Fig. 2), triggering plate bending in the vicinity of the trench (x = 1110 km). The thermal, mechanical and compositional structure of the subduction channel is then free to evolve. To sustain a constant age at trench of the subducting plate, labelled A t , a specific thermal condition is imposed in the top left hand corner of the simulation box. A constant plate age is imposed 710 km away ahead of the trench, by setting on a 400 km wide domain a plate cooling from asthenosphere outcropping at x = 0 to thermal thickening corresponding to an oceanic plate age, A lith , reached at x = 400 km (Fig. 2). As it takes ∼11 Myr for a column of oceanic plate to drift from x = 400 km to the trench, A lith is set to: A lith = A t -11 (in Myr). Crustal tracers are continuously supplied down to the Moho to renew the layer covering the subducting plate in the top left corner. The Moho depth is either ∼7 or ∼8 km depending on resolution mesh (Table 1). Other kinematical, mechanical, and thermal conditions are summarised in Figure 2. To estimate the upper plate stress state, the regional horizontal tectonic force, F res , is computed at x = 1665 km (Fig. 2) from σ xx deviation from lithostatic stresses in an oceanic column of density ρ ref [START_REF] Christensen | An Eulerian technique for thermomechanical modeling[END_REF]:

F res = - zi 0 σ xx dz + zi 0 ρ ref (z)gzdz
, where z i is the compensation depth (278 km) and σ xx is the horizontal normal component of stress tensor. F res corresponds to the force applied to keep constant convergence velocity.

Rheology

A non-Newtonian viscous rheology and a pseudo-brittle rheology are combined. Strength , for both oceanic crust and mantle rocks. The non-Newtonian ductile viscosity, ν d , depending on depth (z), temperature (T ), and second invariant of the deviatoric strain rate tensor ( ε), writes as:

ν d = A 0 exp E a (C) + V a ρgz nRT ε1/n-1 (1)
where A 0 is a pre-exponential factor, E a is the activation energy depending on the lithology (C), V a is the activation volume, P is the lithostatic pressure, n is a exponent higher than 1, ρ is the mantle density, g is the gravity acceleration, and R is the ideal gas constant. In the brittle domain, the yield stress, τ y , writes as :

τ y = τ 0 + γ(C)ρgz (2)
where τ 0 is the frictional cohesive strength at the surface, γ is a parameter related to the friction coefficient [START_REF] Turcotte | Geodynamics: Applications of continuum physics to geological problems, 2nd Edition[END_REF] depending on lithology. The brittle strain rate, εb , is defined as a function of stress, τ : εb = (τ /τ y ) np εref , where n p is a stress exponent higher than 10 and εref is reference strain rate. Using the yield stress and the brittle deformation rate, an equivalent viscosity in the brittle domain, ν b , is defined as:

ν b = τ y   ε 1 np -1 b ε 1 np ref   (3)
The effective viscosity, ν ef f , is computed assuming that the total strain rate is the sum of brittle and ductile strain rates:

ν ef f = 1/(ν -1 d + ν -1 b ).
The crust is supposed to be much weaker than mantle both in the brittle and ductile realms, to promote strain localization within the interplate crustal layer (Fig. 2). The low crustal friction coefficient (γ c = 0.05) belongs to the classical range used to simulate realistic subduction dynamics (e.g., [START_REF] Tackley | Self-consistent generation of tectonics plates in timedependent, three-dimensional mantle convection simulation. 1. Pseudoplastic yielding[END_REF][START_REF] Bird | An updated digital model of plate boundaries[END_REF][START_REF] Gerya | Why is terrestrial subduction one-sided?[END_REF][START_REF] Van Heck | Planforms of self-consistently generated plates in 3d spherical geometry[END_REF]. Other constants and parameters are listed in Table 2.

Numerical resolution and investigated parameters

The influence of mesh refinement is investigated by testing three numerical meshes, with an increasing refinement of the smallest mesh. The refined domain is limited to the rectangular area that encompasses the subduction plane and the mantle wedge tip (Table 1) where gradients of temperature, strain and stress are maximum. For the lowest resolution ("LR"), minimum node spacings are 2.85 and 2.27 km in the horizontal (x)-and vertical (z)-directions. Both are decreased to 1.26 km in the most refined mesh ("HR"). A third mesh of intermediate resolution (labelled "MHR") is also tested. In all cases, the tracer density, uniform all over the simulation box, respects a minimum number of 7 markers in the smallest cell, which offers a sufficient accuracy in resolving the lithological difference [START_REF] Arcay | Numerical simulations of subduction zones: Effect of slab dehydration on the mantle wedge dynamics[END_REF]. Therefore, the total number of tracers increases from 1.2544×10 6 for grid LR to 5.76×10 6 for grid HR.

The numerical mesh refinement will be first investigated. The influence of the subducting crust properties is tested nex. As mantle characteristics are fixed, I focus on the offset in activation energy between crust and mantle, ∆E a : ∆E a = E m a -E c a , mantle being always stiffer than crust. The impact of ∆E a on the down-dip extent of the subduction channel will be explored in details for two end-member ages at the trench of the subducting lithosphere (20 Myr and 100 Myr). All other parameters in equation 1 are kept identical. The density difference between crust (ρ c ) and mantle (ρ m ), ∆ρ c (=ρ c -ρ m ), will be tested only for a young (20 Myr old) subducting plate, only moderately negatively buoyant. I will show how the pair (∆E a , ∆ρ c ) modify the regime of young plate subduction. As a thick (100 km) upper plate is an usual feature of subduction zones but also an unfavorable condition for subduction initiation, it is imposed in all simulations.

Estimate of the interplate down-dip extent

The down-dip extent of the subduction channel may be considered either as a maximum depth of kinematic decoupling [START_REF] Furukawa | Depth of the decoupling plate interface and thermal structure under arcs[END_REF]Arcay et al., 2007b;[START_REF] Wada | Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization[END_REF][START_REF] Wada | Common depth of slab-mantle decoupling: Reconciling diversity and uniformity of subduction zones[END_REF][START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF], or as a minimum depth of viscous and kinematic coupling between subducting slab and asthenosphere. Since partial kinematic coupling might occur at depth shallower than the asthenosphere (e.g., [START_REF] Wada | Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization[END_REF][START_REF] Syracuse | The global range of subduction zones thermal models[END_REF]van Keken et al., 2011, Fig. 1), it is more convenient to tag the subduction channel down-dip extent as a depth of significant slab-asthenosphere coupling, z coup . To estimate z coup , let us zoom in the interplate structure. At convergence start, as the incipient subducting slab enters the asthenosphere, the corner flow is directed towards the mantle wedge tip before turning downwards, dragged by the subducting slab (Fig. 3a,2.35 Myr). The following cooling of the incoming plate rises thermal gradients at the mantle wedge tip, which enhances the corner flow and slightly shallows z coup (Fig. 3a,4.4 Myr). The thermal structure and its associated z coup are stable (Fig. 3a,25.3 Myr). Along the subduction channel, stress (τ ) and strain ( ε), and thence rate of energy dissipation (τ × ε), are very high up to the location where high strain rate jumps into the asthenosphere, where viscosity is low (Fig. 3d). The interplate geometry can hence be tracked by a x -z profile joining maximum dissipation rates modelled at depth shallower than ∼120 km [START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF]. Along this profile, the average strain rate, ε, is close to 2.6×10 -13 s -1 (Fig. 3d1), i.e., ∼ v c /h c , where v c is convergence velocity (6.5 cm/yr) and h c is the weak layer thickness. The interplate geotherm increases roughly linearly with depth until a high temperature increase that marks the proximity to asthenosphere (Fig. 3d2). Therefore, z coup is defined by the depth along the interplate profile where a maximum thermal gradient and a strong strain rate decrease occur. Note that the distribution of high deformation rates, here obtained with a free slip-non deforming top surface and a weak crustal layer, reproduces well the strain rate pattern obtained in convection models including a free surface and a weak crust that simulate self-consistent singled-sided subduction [START_REF] Schmeling | A benchmark comparison of spontaneous subduction models -Tonwards a free surface[END_REF] [START_REF] Quinquis | The role of boundary conditions in numerical models of subduction dynamics[END_REF][START_REF] Crameri | Parameters controlling dynamically selfconsistent plate tectonics and single-sided subduction in global models of mantle convection[END_REF].

Influence of the numerical mesh

A preliminary set of experiments with a ∼100 Myr old, ∼100 km thick subducting lithosphere is performed to study how the numerical resolution affects the maximum interplate decoupling depth, z coup (simulations S1-10a to d, Table 3). The physical properties of the crust are kept constant. Using a coarse mesh (labelled LR, x-and z-spacings of 2.8 and 2.3 km, Table 1) yields a shallow z coup , that stabilizes at the location where the crust content at computational nodes is lower than 99.9%. To understand this result, note that the weak layer slightly thins, probably as a result of slab bending and/or internal shearing. Consequently, the weak layer is discretized by less than 4 nodes by the Eulerian mesh at ∼80 km depth (Figure 4a,top). Close to the crust surface the average crust content at nodes drops from 100% to less than 99.1%. A "contamination" by mantle tracers lower than 1% is sufficient to notably increase the strength averaged at nodes, which triggers the strain rate jump towards asthenosphere at low depth (80 km,Figure 4a,top). With a more refined mesh (HRlim, xand z-spacings of 1.26 km), the weak crust is discretized by more than 4 Eulerian nodes, among which 3 nodes at least display crust contents close to 100%: this enables the computation of an effective low strength at "crustal" nodes and deepens z coup to the level where the Eulerian mesh is not refined anymore (102 km depth, Figure 4a, bottom). By increasing the downwards extent of the refined area, z coup becomes independent of the numerical mesh (114 km, Figure 4b). A minimum number of 3 nodes to discretise the weak layer core, where crust content equals 100%, appears to be the condition to simulate a mesh-independent coupling depth, as also observed by [START_REF] Crameri | A free plate surface and weak oceanic crust produce single-sided subduction on Earth[END_REF]. Refining the mesh increases node number but also the number of compositional tracers, to respect a minimum of 7 particles in the smallest mesh (tracer number multipleed by ∼5, Table 1). The increase in node number and tracer density between grids LR and HR considerably rises computation time. To find a balance between discretization quality and computation time, a third numerical mesh is tested, with intermediate node spacings in the subduction plane area (MHR mesh, node spacing of ∼1.7 km). To improve crust discretization, the crustal layer is also slightly thickened (from 6.9 to 8.3 km, Table 1, Figure 4c). The z coup depth obtained with mesh MHR (109 km) is not exactly the one simulated with mesh HR (114 km). However, ∼85% of the difference between high and low resolution meshes is accounted for with the intermediate mesh, while mesh MHR allows for saving noticeably computation time. Mesh MHR is used in all simulations presented in the following. Numerical accuracy will be further discussed in section 6.1.

Modelling results

Case 1: Subduction of a cold lithosphere

In simulations S1-2 to S6-15 (Table 3), only the offset in activation energy between crust and mantle, ∆E a , is varied, from 20 to 155 kJ/mol (E c a varying from 445 to 310 kJ/mol). This large range of crust viscous strength enables to test different hypotheses on the lithology forming the subduction channel (sediments, basalt, gabbro, mantle, altered/hydrated phases, composite material). The lowest activation energy (310 kJ/mol, ∆E a =155 kJ/mol) mimics a ductile strength encompassed between the ones of a dry quartzite [START_REF] Ranalli | Rheological stratification of the lithosphere[END_REF] and of a wet diorite (Ranalli, 1995, Fig. 5a). The strongest rheology (E c a =445 kJ/mol, ∆E a =20 kJ/mol) simulates the strength of a dry mafic granulite [START_REF] Wilks | Rheology of some continental lower crustal rocks[END_REF].

Influence of interplate strength on the corner flow structure

Softening the interplate channel (increasing ∆E a ) increases z coup (Figure 5b). The largest z coup , 163 km depth, is modelled for ∆E a =135 kJ/mol. Conversely, the stiffest subducting crust (∆E a =20 kJ/mol) shallows z coup to 60 km depth. This ∆E a -z coup trend was already observed by [START_REF] Wada | Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization[END_REF] and described in [START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF], but the values of z coup obtained here are significantly different thanks to the improved resolution. Surprisingly, for the weakest crust (∆E a =155 kJ/mol), the jump of localized strain from the interplate to the asthenosphere does not occur in the numerically refined domain (Fig. 6a2): The crust is so weak that it remains softer than the asthenosphere, which maintains strain localization within the subducting crust. z coup is shifted at the boundary between fine and coarse meshes, at 250 km depth. As the warm corner flow is offset downwards, the fore-arc lithosphere significantly thickens by cooling. The major part of the cold mantle wedge is motionless. For 135< ∆E a <155 kJ/mol, an in-between behaviour is modelled. Near the base of the overriding plate (∼100 km depth), convergence is partly transmitted to the fore-arc base (velocities < 0.5 cm/yr, Fig. 6a1). I label z partial-coup the shallowest depth along the subduction channel where partial kinematic coupling is modelled (Fig. 6a1). It equals 110±10 km in simulation S1-14 (∆E a =145 kJ/mol). Partial kinematic coupling between slab and asthenosphere progressively increases downwards. z coup is defined instead by the coupling depth where rocks overlying the downgoing plate are asthenospheric (hot temperature) and are dragged with velocities close to subduction speed (> 3 cm/yr). According to this definition, z coup for ∆E a = 145 kJ/mol varies through time but globally deepens, up to 200 km depth (Figure 5b). The z coup dependence on time is likely to be a numerical effect, related to transient thinnings of the weak layer that would possibly disappear by increasing resolution. Nonetheless, the layer of partial kinematic decoupling forms on the whole at the upper plate-asthenosphere boundary, independently of numerical resolution. Note that the ∆E a range resulting in partial kinematic coupling is quite narrow (135 < ∆E a < 155 kJ/mol). For a bit younger subducting plate (A t =60 Myr, simulations S6-2 to S6-15), the dependence of the interplate structure on crust strength is similar to the one simulated for a 100 Myr old lithosphere (Figure 5b). The main differences are subduction failure if the oceanic crust is too resistant (∆E a = 20 kJ/mol) and a general slight z coup deepening compared to A t =100 Myr. The striking cooling of the fore-arc lithosphere is modelled when ∆E a ≥ 145 kJ/mol, and partial kinematic coupling occurs at the base of the fore-arc since ∆E a = 135 kJ/mol (Figure 5b). The numerically-induced variations in z coup are however much wider than simulated when A t =100 Myr. To summarise, the fore-arc base heating or cooling depends on the corner flow structure, which depends on ∆E a . The downwelling stream of the corner flow is well developed whatever crust strength, since the asthenosphere sticking to the slab ("viscous blanket", Kincaid and Sacks, 1997, Fig. 3b) does not depend on ∆E a . On the contrary, the horizontal part of the corner flow is fast and focussed at the upper plate base only if ∆E a is rather low (< 135 kJ/mol, Fig. 3a), and provides in this case a significant advected heat at the fore-arc base that thins it. The horizontal mantle flow gets more and more diffuse as ∆E a is risen, which restrains fore-arc basal heating (compare Fig. 6a to Fig. 3a). When ∆E a = 155 kJ/mol, the heat advected towards mantle wedge tip is so reduced that the fore-arc conductive cooling is unbalanced.

Effect of a very weak subduction channel on partial melting in the mantle

wedge To explore the consequences of an extreme cooling of the fore-arc lithosphere when ∆E a ≥145 kJ/mol, water transfers associated with slab dehydration and overlying rock hydration are computed in simulations S1-14 and S1-15 (Fig. 6a3). The subducting oceanic crust and the underlying layer of serpentinised mantle (7 km thick) are both assumed to be water-saturated before subduction (containing respectively 2.93 and 6.5% in weight of water). Phase transitions controlling water transfers are predicted by using pressure-temperature paths of lithological tracers, based on accurate phase diagrams determined for a water-saturated peridotite [START_REF] Schmidt | Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation[END_REF] and a gabbroic composition [START_REF] Bousquet | Kinematic, thermal and petrological model of the Central Alps: Lepontine metamorphism in the upper crust and eclogitisation of the lower crust[END_REF]. The model of water expelling, upwards migration, and absorption by the nearest under-saturated marker above the dehydrating element, extensively described and tested in previous papers [START_REF] Arcay | Numerical simulations of subduction zones: Effect of slab dehydration on the mantle wedge dynamics[END_REF][START_REF] Arcay | Overriding plate thinning in subduction zones: Localized convection induced by slab dehydration[END_REF](Arcay et al., , 2007b[START_REF] Arcay | Back-arc strain in subduction zones: Statistical observations vs. numerical modelling[END_REF] will not be detailed here. Water transfers are passive and do not affect mantle and crust properties. The key result of the experiment S1-15 is that the mantle wedge hydrated by slab dehydration is too cold to allow for the triggering of mantle melting in wet conditions. Let us assume that the lowest temperature of wet mantle solidus is 980 o C [START_REF] Green | Experimental melting studies on a model mantle composition at high pressure under water-undersaturated conditions[END_REF][START_REF] Millholen | Melting interval of peridotite with 5.7 per cent water to 30 kilobars[END_REF][START_REF] Schmidt | Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation[END_REF]. The progressive cooling of the fore-arc lithosphere shifts away from the trench this critical isotherm (Fig. 6a3): at times longer than 16.4 Myr, all the water released by slab dehydration saturates a mantle wedge whose temperature is well below 980 o C, as above one slab dehydration front a column of chlorite forms (T 850 o C, water content of 2.1 %wt, Fig. 6a3,inset). Consequently, the model predicts that the volcanic arc activity of the subduction zone should stop despite a significant fluid flow ex-pelled from the slab, which is the eventual outcome of a very weak subduction channel. Regarding anhydrous mantle melting, it is all the more impeded as it required temperature hotter than wet partial melting. When ∆E a =145 kJ/mol (simulation S1-14), the fore-arc cooling is not as excessive as in the preceeding case: Partial melting within the hydrated mantle wedge is allowed in a limited hot region. For A t =60 Myr (simulation S6-13), partial kinematic coupling enhances fore-arc cooling through time. The mantle located ontop of the slab dehydration front the farthest from the trench (chlorite out) is predicted to be too cold to enable mantle melting (T ≤ 980 o C) when convergence duration exceeds 20 Myr. The use of the wet solidus determined by [START_REF] Green | Experimental melting studies on a model mantle composition at high pressure under water-undersaturated conditions[END_REF] and [START_REF] Millholen | Melting interval of peridotite with 5.7 per cent water to 30 kilobars[END_REF] to assess mantle melting occurence in water-saturated condition will be discussed in section 6.

Case 2: Subduction of a young and warm lithosphere

In this section, experiments simulate the subduction of a young incoming lithosphere, ∼40 km thick and aged of 20 Myr at trench. For a hot downgoing plate, the efficiency of mechanical decoupling across the subduction interface may be a function not only of the subduction channel strength (∆E a ) but also of the weak layer density (∆ρ c ). Indeed, the total slab weight favouring subduction is more sensitive to crust density than in the case of a cold incoming plate, whose mantle density is high and lower the relative effect of crust density. The range of investigated ∆E a is enlarged to 220 kJ/mol (E c a =245 kJ/mol), that simulates the weak strength of a granite (Mackwell et al., 1998, Fig. 5a). The density of the oceanic crust varies between 2920 kg.m -3 (∆ρ c = -380 kg.m -3 ), which is close to the density of gabbro (2940 kg.m -3 ) or andesite (2900 kg.m -3 ) in the greenschist facies [START_REF] Bousquet | Kinematic, thermal and petrological model of the Central Alps: Lepontine metamorphism in the upper crust and eclogitisation of the lower crust[END_REF], and 3550 kg.m -3 (∆ρ c = +250 kg.m -3 ), corresponding to the density of an eclogitized gabbro [START_REF] Bousquet | Kinematic, thermal and petrological model of the Central Alps: Lepontine metamorphism in the upper crust and eclogitisation of the lower crust[END_REF]. The mantle density is kept constant (3300 kg.m -3 at the surface).

Subduction regime as a funtion of ∆E a and ∆ρ c

Figure 7 sums up the broad variety of subduction patterns modelled for A t =20 Myr depending on the pair (∆E a , ∆ρ c ) (simulations S2-2a to S2-22g, Table 3). Different behaviours may be modelled simultaneously, yielding overlappings domains. Domain 1 may be defined where subduction aborts rapidly (<10 Myr): the subduction plane is blocked and convergence localises ahead of the trench (white domain in Fig. 7, inset 1). For a stiff crust (∆E a ≤ 40 kJ/mol), subduction is possible only if crust is denser than the mantle (∆ρ c >0). The softer the crust is, the wider the range of crust density allowing subduction is (yellow domain). However, near the outline of the subduction domain, stresses within the upper plate are highly compressive (light brown domain in Fig. 7, inset 4): subduction would probably stop if convergence was driven by a constant force less than 10 13 N/m (computed as explained in section 3.1) instead of constant velocity. When the crust is denser than the mantle, early slab breakoff occurs if crust is too stiff, followed by subduction obstruction (∆E a 40 kJ/mol, left green domain, inset 6). At slightly higher ∆E a and positive ∆ρ c , the incipient subduction zone is unstable, and the slab finally detachs from the surface (purple domain, inset 7). On the contrary, for a weak and dense crust, the early slab break-off does not repeat and subduction is stable (∆E a 105 kJ/mol, right green domain, inset 5). The initial slab drip could be related to the weakness of the incipient plate [START_REF] Andrews | Rheologic controls on the dynamics of slab detachment[END_REF], and/or to the low activation volume [START_REF] Crameri | A free plate surface and weak oceanic crust produce single-sided subduction on Earth[END_REF]. When the crust is weak and buoyant, small crust plumes initiate near the upper lithosphere-asthenosphere boundary (∼110±15 km depth, if the crustal geotherm at this depth is hot enough (T ∼ 1000 o C, Fig. 6b2; red domain in Fig. 7 inset 2). Two levels of crust plume formation are modelled in a few simulations, the deeper one occuring at ∼130 km depth (section 5.2.2). Crust plumes have been modeled by [START_REF] Vogt | Numerical modeling of geochemical variations caused by crustal relamination[END_REF] and [START_REF] Maunder | Relamination of mafic subducting crust throughout earth's history[END_REF] in the same depth range for similar densities and rheological parameters, but for a much thinner upper plate (∼50 km) and for slower convergence. In the present study, crust diapirism may be promoted by shear heating in the subduction channel. In the Archean, the underplating of crust plumes along the overlapping plate Moho ("relamination") has been suggested to explain the lower continental crust differentiation, as crust plumes may include partially molten crust [START_REF] Hacker | Differentiation of the continental crust by relamination[END_REF]. Crust plumes do not form for an extremely weak and buoyant crust (∆E a = 220 kJ/mol and ∆ρ c = -380 kg.m -3 , simulation S2-22a). The subducting crust thickens the subduction channel by accumulating at depth shallower than ∼100 km and locally convect (velocity of ∼5 cm/yr, blue domain in Fig. 7, inset 3). The subduction plane thickening, modeled by [START_REF] Maunder | Relamination of mafic subducting crust throughout earth's history[END_REF] for a similar crustal rheology, was deemed to be a plate tectonics regime in the Archean [START_REF] Hoffman | Archean oceanic flake tectonics[END_REF].

Partial kinematic coupling and thermal state of the fore-arc lithosphere

As the formation of crustal plumes triggers partial kinematic coupling at the extremity of the subduction channel, the experiments showing crustal diapirism are detailed in this section, illustrated by simulation S2-17b. When instabilities of buoyant crust initiate at ∼100 km depth, z coup is shifted to 137 km depth while partial kinematic coupling occurs between the plume root and z coup (Fig. 8a, 9 Myr). Soon a secondary branch of asthenospheric corner flow forms reaching the root of the main crust plume (∼90 km depth, 14.2 Myr). This rises z coup to z partial-coup , leading to a classical velocity field with a fast corner flow (>3 cm/yr, Fig. 8a, 15.2 Myr), before z coup deepens again to ∼190 km at 18.9 Myr. The crust flow within the main plume (clock-wise current) cools the fore-arc. The crust spinning sense is reversed when z coup rises up to the crust plume head, resulting in a fore-arc warming (Fig. 8, 14.2 Myr). The sequence of forearc cooling induced by crust upwelling, followed by a fast warming related to asthenospheric rising, occurs repeatedly. Crust instability thus basically entails (1) a high time-dependence of the fore-arc structure, and (2) an intermittent but intense cooling of the fore-arc domain. Outside the domain of crust instability, partial kinematic coupling appears when crust buoyancy exceeds -140 kg.m -3 at ∼114 km depth. Figure 8b depicts the typical evolution simulated in simulation S2-13g: z coup deepens with time (from 175 km to ∼200 km depth between 8.1 and 10.2 Myr), before being suddenly shallowed when the asthenosphere rises up to z partial-coup (113 km,12,8 Myr). The z coup deepening followed by quick shallowing is once again a repeating process (Fig. 8).

z coup dependence on ∆E a , ∆ρ c , and time

Figure 9 displays z coup variations with crust strength and buoyancy after ∼10 Myr of convergence (panel a). z coup is mainly constant through time in the absence of crust plume, but becomes strongly time-dependent otherwise. The transient state resulting from fast z coup shallowings up to the depth of z partial-coup is captured in panel b, where the shallowest z coup value simulated between 10 and 13.5 Myr of convergence is mapped. First, the large range of z coup modelled for A t =20 Myr varies from 50 km to 250 km depth. This maximum depth modelled at 13.5 Myr in simulation S2-17c (∆E a = 175 kJ/mol; ∆ρ c = -240 kg.m -3 ) is controlled by the mesh coarsening. Second, in the domain where z coup is stable (for ∆ρ c -120 kg.m -3 ), the coupling depth deepens when ∆E a increases, as already described for A t = 100 Myr (section 5.1.1). At constant crust strength, z coup slightly deepens with increasing crust buoyancy, possibly because of a decrease in slab dip angle (compare modelled slab dips in insets 2, 4 and 5 in Figure 7), though the relationship between slab dip angle and crust buoyancy may not be so straightforward (compare insets 5 and 3). Third, when z coup is time-dependent, crust plumes and partial kinematic coupling strongly increase z coup (> 150 km), this deepening being partly counterbalanced by punctual upwards z coup jumps that temporarily reverse the z coup -∆E a trend observed for ∆ρ c > -120 kg.m -3 (solid black squares in Fig. 9b).

Wet mantle melting as a function of z coup evolution

Water transfers are computed in 10 experiments of young plate subduction to evaluate the likelihood of wet mantle melting (purple symbols in Fig. 9a). As mantle melting is not directly computed, I focus on the water flux that should trigger melting. To trigger melting at the slab surface, the water flux expelled by slab dehydration must be not only significant but also absorbed at temperature hotter than 980 o C (lowest wet solidus for peridotite, section 5.1.2). In the present modelling this water is first dissolved in the asthenosphere and fore-arc base (water content at saturation < 0.05 wt%, [START_REF] Arcay | Numerical simulations of subduction zones: Effect of slab dehydration on the mantle wedge dynamics[END_REF]. If there is some left, the remaining water hydrates the cold overlapping lithosphere where nominally hydrated phases form (amphibole, talc, serpentine), increasing water contents above 0.5 wt%. Water contents in the upper plate, that here acts as a "sponge", are thus tracked to estimate the water amount available at slab surface at T > 980 o C that would normally yield mantle melting (Fig. 10). Three melting regimes are modelled. If fore-arc cooling is limited (z coup 140 km), one part of dehydration water should always be absorbed at T > 980 o C in amounts greater than 0.5 wt%, which might be sufficient for sustaining the volcanic front activity even if the corresponding metasomatized domain is narrow (Fig. 10a, solid triangles in Fig. 9). If fore-arc cooling is intense (z coup always ≥ 160 km), almost all the released water is absorbed at temperature below wet solidus (Fig. 10c, crosses in Fig. 9). Therefore, the volcanic front activity would probably vanish. An intermediate regime occurs if mantle wedge cooling is disturbed by sudden z coup shallowings. Compared to the previous situation, transient forearc warmings extend mantle melting for a while (Fig. 10b, empty triangles in Fig. 9). Magma supply and volcanic front activity might be intermittent in this case. In conclusion, the mantle wedge cooling resulting from the subduction channel weakness significantly restricts wet mantle melting for A t = 20 Myr, as already described for older plate subduction.

Synthesis

5.3.1. Decoupling along the subduction interface as a function of slab age, crust buoyancy, and interplate rheology Figure 11 schematizes how crust strength and slab weight regulate kinematic decoupling between converging plates in this modelling. As already mentioned, the weight of the downgoing slab depends here on crust density and on lithosphere age at trench [START_REF] Turcotte | Geodynamics: Applications of continuum physics to geological problems, 2nd Edition[END_REF]. If A t ≥60 Myr, the average slab density exceeds the asthenospheric one whatever crust buoyancy. The resulting slab weight excess promotes high dip angle (> 30 o ) when the incipient slab enters the asthenosphere. This activates the corner flow that eases kinematic decoupling at the interplate down-dip extent and, finally, subduction. This process is only slightly sensitive to the subduction channel strength (purple domain in Fig. 11). Nonetheless, the latter controls the structure of the horizontal branch of the corner flow: fast, focussed and warming the fore-arc base when ∆E a is low (red area in Fig. 12), and inversely slow, diffuse and advecting a low heat flux at the mantle wedge tip when ∆E a is high (blue area). On the contrary, if A t =20 Myr, as mantle slab density is minimum, crust buoyancy has a significant influence on the total slab weight excess. If ∆ρ c = -380 kg.m -3 , the low dip angle of the incipient slab hampers corner flow activation. Subduction initiates only if the interplate decoupling is very efficient up to asthenospheric levels, i.e., if crust is quite weak (yellow domain in Fig. 11). If the subducted crust is not so buoyant (∆ρ c ≥0), the increase in slab dip angle helps the corner flow triggering and lowers the minimum ∆E a enabling subduction. 5.3.2. Fore-arc structure as a function of subducting plate age and crust weakness Figure 12 sums up the major characteristics of the fore-arc structure modelled as a function of incoming plate age, A t , and crust weakness, ∆E a , for a crust buoyancy set to -380 kg.m -3 . To begin with, along the boundary delimiting subduction, the depth of slab-asthenosphere coupling is the shallowest that can be modelled for a given A t , since ∆E a is the lowest (Fig. 5b). Therefore, the fore-arc lithosphere base is predicted to be there the warmest (red area). This is nonetheless not true for a 20 Myr old incoming lithosphere subducting only for very weak crust, which leads to fore-arc cooling and cold crust diapirism (green area). Fore-arc cooling by z coup deepening occurs if ∆E a ≥ 155 kJ/mol whatever subducting plate age, and can be intense enough to prevent mantle melting in wet conditions (blue area). Lastly, partial kinematic coupling at fore-arc base is modelled in a narrow range of crust strength (135 < ∆E a < 155 kJ/mol if A t = 100 Myr), and moderately cools the fore-arc. The depth where partial kinematic coupling starts, z partial-coup , is mainly a function of the upper plate thickness. As partial kinematic coupling is systematically associated with Rayleigh-Taylor crustal instabilities when A t =20 Myr, the domain of partial kinematic coupling for A t ≤ 60 Myr is assumed to merge with the domain of cold crust plumes. To conclude, no universal subduction plane strength is able to simulate alone realistic subduction, whatever subducting plate age, if the crust is significantly buoyant and the overlapping plate thick (∼100 km).

Discussion and concluding remarks

Study limitations

Two main shortcomings can be underlined: (1) the numerical scheme used to interpolate strength from lithological markers'properties, and (2) modelling the subduction channel by a rheological contrast constant with depth. The viscosity at Eulerian nodes can be computed from markers' values (1) using an arithmetic averaging, assuming deformation at constant strain rate (such as pure shear parallel to compositional layering, [START_REF] Schmeling | A benchmark comparison of spontaneous subduction models -Tonwards a free surface[END_REF]. On the contrary, a harmonic averaging simulates simple shear (at constant stress) parallel to the channel surface. While arithmetic averaging yields a numerical strengthening, the harmonic one reproduces well a thin channel low viscosity, and is more adapted to simulate a subduction channel [START_REF] Schmeling | A benchmark comparison of spontaneous subduction models -Tonwards a free surface[END_REF][START_REF] Deubelbeiss | Comparison of Eulerian and Lagrangian numerical techniques for the Stokes equations in the presence of strongly varying viscosity[END_REF]. As both kinds of shear are expected in a subducting slab, the averaging scheme should change depending on deformation pattern, which is not easy to handle numerically [START_REF] Schmeling | A benchmark comparison of spontaneous subduction models -Tonwards a free surface[END_REF]. Geometric averaging leads to a in-between averaged strength. In the present study, brittle and ductile rheologies are combined through harmonic averaging (section 3.2) while viscosity at nodes is estimated from compositional markers by an arithmetic averaging. The subduction channel strengthening simulated for a coarse mesh (section 4) directly results from arithmetic averaging. To better model viscosity layering associated with lithological boundaries, resolution may be improved by either (i) increasing the number of Eulerian nodes or thickening the low viscosity layer when arithmetic averaging is used (as in this study), (ii) decreasing the interpolation volume around Eulerian nodes [START_REF] Duretz | Discretization errors and free surface stabilization in the finite difference and marker-in-cell method for applied geodynamics: A numerical study[END_REF], (iii) limiting irregularities in marker distribution by implementing a conservative velocity interpolation scheme [START_REF] Wang | Advantages of a conservative velocity interpolation (CVI) scheme for particle-in-cell methods with application in geodynamic modeling[END_REF], or (iv) modeling compositional interface during marker advection, through for instance the particle level set method [START_REF] Samuel | Modeling advection in geophysical flows with particle level sets[END_REF].

Imposing a ductile strength contrast ∆E a constant with depth along the subduction channel (2) supposes that lithological variations are small. It is more likely that the medium composing the interplate changes with depth, from, possibly, sediments at shallow depth, upper then lower crusts, to more mafic composition at great depth. This lithological sequence could be simply modelled by a strength increase along the subduction channel, i.e., by a ∆E a decrease with depth (Fig. 5a). As low ∆E a yields z coup shallowing (Fig. 5b and9), a ∆E a lowering with depth should shallow z coup . If softening minerals locally form, such as serpentine at high depth [START_REF] Duarte | How weak is the subduction zone interface[END_REF], a localized ∆E a increase might favor z coup deepening, but this would have to remain limited, otherwise volcanic arc extinction by mantle wedge cooling might occur (Fig. 12).

Partial kinematic coupling

Partial kinematic coupling next to the mantle wedge tip would explain the extent of the "cold fore-arc nose" detected in seismic tomographies (e.g., [START_REF] Zhao | Tomographic imaging of P -and S-wave velocity structure beneath Northeastern Japan[END_REF][START_REF] Zhao | Seismological structure of subduction zones and its implications for arc magmatism and dynamics[END_REF][START_REF] Stachnik | Seismic attenuation and mantle wedge temperatures in the Alaska subduction zone[END_REF]Abers et al., 2006, Fig. 3). Partial kinematic coupling would result from ductile strain localized within weak materials, such as sediments or hydrated phases (serpentine, talc, [START_REF] Peacock | Hydrous minerals in the mantle wedge and the maximum depth of subduction thrust earthquakes[END_REF], or brittle deformation occurring deeper than the brittle-ductile transition. In thermo-kinematic models, partial kinematic coupling is modelled by imposing a given percent of kinematic coupling over a specified length at the subduction plane extremity. Coupling abruptly passes from partial to full to mimic the effect of a dehydration reaction on viscous coupling [START_REF] Kneller | B-type olivine fabric in the mantle wedge: Insights from high-resolution non-Newtonian subduction zone models[END_REF][START_REF] Kneller | Stress, strain, and B-type olivine fabric in the fore-arc mantle: Sensitivity tests using high-resolution steady-state subduction zone models[END_REF][START_REF] Syracuse | The global range of subduction zones thermal models[END_REF]. A low coupling (2-3%) between 40 and 70 km depth cools the fore-arc while sustaining a sharp sub-vertical boundary between low-temperature fore-arc and hot asthenospheric wedge, consistent with seismic cross-sections and heat flow profiles. An increased coupling (∼10%) leads to a cooling of the subduction plane down-dip extent while a "hot tongue" of asthenosphere rises up to the location where full decoupling between converging plates is imposed (∼40 km depth). This odd thermal structure is actually not consistent with seismic tomographies. The present study shows that partial kinematic coupling may form according to the incoming plate age and the subduction plane strength (135 -145 kJ/mol if A t =100 Myr, enlarged to ∆E a > 120 kJ/mol if A t =20 Myr, Fig. 7). The transition to full kinematic coupling is then progressive (from 0 to ∼7%) and occurs on a large distance (> 100 km extent), while it is pinpointed and thus very sharp in thermo-kinematic models. The thermal structure modelled at slab surface is then much more smoothed than in thermo-kinematic modellings.

Diapirism-induced fore-arc cooling, numerical resolution, and mantle melting prediction

In this modelling, z coup is mesh-independent if the weak layer core is discretized by more than 3 computational nodes. Otherwise, the artificial strengthening of the subduction channel triggers a hot mantle upwelling to soften the interplate, ending in a numerical shallowing of z coup . When crust plumes are modelled, they reduce the crustal channel thickness. As a consequence, the asthenospheric corner flow frequently rises up to the root of crust instabilities. The sporadic upwards jumps of z coup are thus probably a numerical effect that could be limited using a more refined mesh. The thermal state of the fore-arc would then be even colder and wet mantle melting even more hindered than simulated for A t =20 Myr. This is consistent with the general fore-arc cooling modelled for older incoming plates at high crust strength reduction.

Can nominally hydrated phases be stable at the slab surface?

Two main layers of hydrous phases are assumed to release the water responsible for partial melting in the mantle wedge. The deeper one is located within the subducting slab, inside the oceanic crust and/or the underlying serpentinized mantle. The second one would be the viscous blanket covering the downgoing slab (Fig. 3b), hydrated by the underlying dehydrating layers. Amphibole and serpentine have been assumed to be stable at the slab surface [START_REF] Iwamori | Numerical modeling of thermal structure, circulation of H 2 O, and magmatism-metamorphism in subduction zones: Implications for evolution of arcs[END_REF][START_REF] Wada | Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization[END_REF][START_REF] Hirauchi | Rheological contrast between serpentine species and implications for slab-mantle wedge decoupling[END_REF], but chlorite would the hydrated phase stable at the hottest temperature at high depth (∼ 860 o C at 60-100 km depth, [START_REF] Schmidt | Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation[END_REF][START_REF] Grove | The influence of H 2 O on mantle wedge melting[END_REF][START_REF] Grove | Kinematic variables and water transport control the formation and location of arc volcanoes[END_REF][START_REF] Till | The beginnings of hydrous mantle wedge melting[END_REF]. Chlorite destabilization would provide the water triggering melting at the bottom of the mantle wedge [START_REF] Till | The beginnings of hydrous mantle wedge melting[END_REF], see section 6.5). Still, chlorite stability was assessed based on thermal models performed with a temperaturedependent Newtonian viscosity. If the stress-dependence of olivine rheology is included, the increase in thermal gradient at the slab surface [START_REF] Van Keken | High-resolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle[END_REF][START_REF] Syracuse | The global range of subduction zones thermal models[END_REF] prevents the formation of stable chlorite. Similarly, the simulations performed in this study predict that chlorite cannot form at the base of the mantle wedge if ∆E a <135 kJ/mol (Fig. 13). Chlorite formation at high depth is predicted when the simulated subduction zones are the coldest, notably when z coup is extremely deep (> 160 km depth, Fig. 6a3 and10a,b). However, since hydrous minerals may be softening materials, their formation would lower the asthenosphere viscous drag at the slab surface, emphasizing wedge cooling and melting hindrance. Finally, as the crust mainly dehydrates before reaching the asthenosphere (Fig. 6a3), this study suggests that the deserpentinization of the subducting mantle may be the major source of the water triggering melting in the mantle wedge.

Solidus assumed to assess wet mantle melting

When the subduction channel is particularly weak (∆E a ≥135 kJ/mol, or a bit more depending on ∆ρ c if A t =20 Myr), the water released by slab dehydration is predicted to be absorbed in the mantle wedge at temperature below the wet solidus, switching off the volcanic front activity. The model assumes the minimum temperature of the wet solidus to be 980 o C, as estimated by [START_REF] Kushiro | Melting of a peridotite nodule at high pressures and high water pressurs[END_REF] and [START_REF] Green | Water and its influence on on the lithosphere-asthenosphere boundary[END_REF]. However, other experiments performed on longer durations [START_REF] Mysen | Melting of a hydrous mantle: I. Phase relations of natural peridotite at high-pressures and temperatures with controlled activities of water, carbon dioxide, and hydrogen[END_REF][START_REF] Grove | The influence of H 2 O on mantle wedge melting[END_REF][START_REF] Till | The beginnings of hydrous mantle wedge melting[END_REF], and/or based on peridotite samples less depleted in iron [START_REF] Grove | The role of H 2 O in subduction zone magmatism[END_REF] suggest a lower wet mantle solidus, related to the chlorite stability field: 810±10 o C between ∼2 and 4 GPa. Therefore, the inhibition of mantle melting induced by an extreme z coup deepening might not occur if the relatively cold wet solidus of [START_REF] Till | The beginnings of hydrous mantle wedge melting[END_REF] was taken into account. Nonetheless, the strong fore-arc cooling (modelled for instance in simulation S1-15, Fig. 6a3) is artificially limited by the numerical mesh coarsening imposed at 250 km depth. The fore-arc would thus be even more cooled if mesh coarsening was shifted to the base of the simulation box. Finally, I assume that, at least on long-time scales (>30 Myr), the fore-arc conductive cooling would prevent mantle partial melting whatever the considered wet solidus.

Slab surface temperature

To discuss the range of realistic subduction channel strengths, let us consider the thermal state modelled along the subducting slab surface, defined as the top of the downgoing oceanic crust, for A t =100 Myr. (Fig 13d). Slab surface temperature (SST) are compared with pressure-temperature (P -T ) conditions estimated from geothermometry derived from H 2 O/Ce measurements on melt inclusions [START_REF] Plank | Emerging geothermometers for estimating slab surface temperatures[END_REF][START_REF] Cooper | Global variations in H 2 O/Ce: 1. Slab surface temperatures beneath volcanic arcs[END_REF]. Data are compiled from very different subduction zones, such as slow subduction of young plates (e.g., Cascadia) or fast subduction of old plates (NE Japan or Tonga). They cannot not be directly compared with modelled SST, but rather considered as upper and lower boundaries. Modelling an overlapping plate thinner than simulated would also probably shift SST to lower depth. The modelling subduction setup simulates the one prevailaing in the Kamchatka [START_REF] Heuret | Plate motions, slab dynamics and back-arc deformation[END_REF]. Imposing a stiff oceanic crust (∆E a =20 kJ/mol, Fig 13a) yields a slab surface much hotter than observed, while the strong fore-arc cooling associated with a weak suduction channel (∆E a ≥135 kJ/mol) produces too cold SST (Fig 13b andd). Imposing ∆E a =80 kJ/mol reproduces the hotter observed SST, whereas SST obtained with ∆E a =105 kJ/mol correctly simulate the P -T conditions observed in the Kamchatka. We should expect the modelled SST to be colder than estimated along several young and/or slow subducting plates (South Chile, Nicaragua or Lesser Antilles), but it is the inverse that is simulated. This would argue for increasing ∆E a to ∼120 kJ/mol to moderate SST increase. I conclude that the ∆E a range that mimics best worldwide SST estimates from thermometry is 80-120 kJ/mol. Moreover, the shallowest and hottest nose in the mantle wedge may be constrained by conditions of last equilibration with the mantle of near-primary arc magmas (Kelemen et al., 2003, Fig. 13a), considered as the signature of hydrated melts reaching the upper plate Moho. Still, hot conditions at the arc Moho may derive from the heat advected by melt migration, which is here neglected. The hottest mantle wedge is modelled for the stiffest subduction channel, but the hot nose is not as shallow as observed (∆E a =20 kJ/mol, lines with stars in panel 13d). To reproduce observations, an additional ingredient must be included to enhance the corner flow. In simulation S1-8b, first, water transfers associated with slab dehydration/mantle wedge hydration are simulated (section 5.1.2); second, an hydrous strength reduction whose amplitude increases with water content [OH -] is imposed for the hydrated mantle: low reducing factor (increasing from 5 to 32) when water is dissolved in nominally anhydrous minerals (200 ppm ≤ [OH -] ≤0.2 wt%), high and maximum (set to 350) when hydrous minerals form (e.g., amphibole, OH -≥0.5 wt%, Arcay et al., 2005, Fig. 13c). ∆E a is set to 80 kJ/mol. Viscosity is then low in the hydrated mantle wedge (∼10 18 Pa.s at ∼50 km depth, Fig 13c), which shallows z coup to ∼48 km and simulates the shallow and hot mantle nose inferred from near-primary arc magmas (1250 o C at ∼50 km depth, Fig. 13d). The drawback of this water-related weakening is that Rayleigh-Taylor instabilities initiate at the base of the hydrated upper plate (Arcay et al., 2005, Fig 13c), thickening the viscous blanket. Compared to simulation S1-8 without water transfers, SST is cooled by ∼150 o C, mimicking the coldest P -T conditions inferred from H 2 O/Ce measurements. As the convective destabilization of the overriding lithosphere is impeded if the maximum hydrous strength reduction is lower than 20 3/2 90 [START_REF] Arcay | Overriding plate thinning in subduction zones: Localized convection induced by slab dehydration[END_REF], a compromise could be to limit the maximum water-related weakening to 90, instead of 350. The asthenospheric viscosity would then be lowered by a factor varying from 4 to 19. 6.7. Interplate strength favouring subduction viability and reproducing cold forearc nose extent This study shows that there is no universal value of the subduction channel strength able to model realistic subduction whatever subducting plate age (i.e., successful initiation, limited compressive stresses and occurence of wet mantle melting), if crust buoyancy is -380 kg.m -3 (Fig. 12). Since the best ∆E a interval modelling reasonable SST is 80-120 kJ/mol (section 6.6), crust buoyancy should not exceed -80/-200 kg.m -3 to allow young plate subduction (Fig. 7). Nevertheless, the simulation of a hot and shallow mantle wedge consistent with conditions of last equilibration of near-primary arc magmas requires a quite different range of crust strength in anhydrous conditions: ∆E a ≤20 kJ/mol. Moreover, the sharp transition between low to high seismic velocities delineating the cold fore-arc nose is often observed at shallow depth: ∼65 km in NE Japan and Alaska [START_REF] Zhao | Tomographic imaging of P -and S-wave velocity structure beneath Northeastern Japan[END_REF][START_REF] Stachnik | Seismic attenuation and mantle wedge temperatures in the Alaska subduction zone[END_REF][START_REF] Abers | The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow[END_REF]; ∼45 km in Hikurangi [START_REF] Davey | Fore-arc mantle wedge seismicity under northeast New Zealand[END_REF] and in the Cascadia [START_REF] Zhao | Seismological structure of subduction zones and its implications for arc magmatism and dynamics[END_REF] and <40 km in S Japan [START_REF] Miura | Structural characteristics controlling the seismicity of southern Japan Trench fore-arc region, revealed by ocean bottom seismographic data[END_REF]. If seismic anomalies image thermal contrasts, to model such shallow z coup when A t 100 Myr, ∆E a should be set to ∼30-40 kJ/mol in Alaska and NE Japan, and lower than 20 kJ/mol in Hikurangi (Fig. 5b); in the Cascadia and S Japan where A t <25 Myr, ∆E a should be 20-40 kJ/mol and at most set to 20 kJ/mol, respectively (Fig. 9). The constraints from seismic tomography lead to 10 ∆E a <40 kJ/mol, which is consistent with the range inferred from hot mantle wedge-derived melts. Finally, two different ranges of ∆E a may be underlined: 80-120 kJ/mol and 10-40 kJ/mol, that are not compatible although constraints from petro-geochemistry and geophysics were obtained on the same subduction zones. How to infer a reconciling value? Again, an extra process regulating the mantle wedge structure could be invoked, such as a decrease in mantle strength depending on [OH -] (section 6.6). A moderate weakening (∼one order of magnitude) in the hydrated wedge should be sufficient to simulate realistic SSTs and a shallow hot wedge core, keeping for the interplate strength the median value ∆E a =80 kJ/mol (simulation S1-8b, Fig. 13c-d).

6.8. Is the down-dip extent of the subduction plane likely to be identical for worldwide subduction zones? In conclusion, three sets of rheological parameters are pointed out: (1) ∆E a ∼ 80-120 kJ/mol to simulate SST derived from H 2 O/Ce measurements; (2) ∆E a ∼10-40 kJ/mol to model the hot and shallow mantle wedge inferred from seismic tomographies and near-primary arc magmas; and (3) a median crust strength (∆E a ∼80 kJ/mol) combined to a water-related weakening. Only rheological sets (1) and (3) are compatible. In case (2), z coup is strongly controlled by the corner flow. As the corner flow dynamics basically depends on the incoming slab characteristics and subduction parameters (such as convergence rate), if ∆E a = 25±15 kJ/mol then the subduction plane down-dip extent very probaly varies from a subduction zone to another. In case (1), z coup equals more or less the upper plate thickness: a bit shallower if ∆E a = 80 kJ/mol (82 km if A t =100 Myr), or a bit deeper if ∆E a = 120 kJ/mol (135 km, Fig. 5b). Hence, the down-dip extent of the subduction channel may be almost independent from subduction parameters if ∆E a = 120 kJ/mol, and moderately dependent if ∆E a =80 kJ/mol. In case (3), z coup strongly depends on subduction velocity (as shown for ∆E a =110 kJ/mol by [START_REF] Arcay | Back-arc strain in subduction zones: Statistical observations vs. numerical modelling[END_REF], using a mesh too coarse implying crust strengthening: the effective ∆E a might have been ∼80 kJ/mol). Since rheology (3) enables to verify all the listed constraints on the mantle wedge thermal state, it may be the best choice to simulate the subduction channel strength. This implies to limit crust buoyancy to ∆ρ c ∼-80 kg.m -3 to allow young plate subduction (Fig. 7), assuming crust eclogitization at ∼17 kbar [START_REF] Bousquet | Kinematic, thermal and petrological model of the Central Alps: Lepontine metamorphism in the upper crust and eclogitisation of the lower crust[END_REF]. Besides, if ∆E a ∼ 80 kJ/mol, then partial kinematic coupling at the fore-arc base should not occur. Lastly, the rheological properties of the subduction channel estimated in this study imply that the subduction plane down-dip extent varies as a function of the subduction set-up, as suggested by the wide range of slab sub-arc depths (England et al., 2004;[START_REF] Syracuse | Global compilation of variations in slab depth beneath arc volcanoes and implications[END_REF], section 1). This conclusion disagrees with a single value of z coup universally equal to 80 km depth [START_REF] Wada | Common depth of slab-mantle decoupling: Reconciling diversity and uniformity of subduction zones[END_REF].

(National Institute of Universe Science) program "TelluS-SYSTER" (2011, 2015 and 2016). The red area depicts the mantle wedge part hot enough for melting in wet conditions. In the bottom panel, two kinds of kinematic coupling between the subducting slab and the overlying mantle are considered: In case of full coupling, velocities at the slab surface are close to the subduction speed as soon as the down-dip extent of the subduction channel is crossed, defining zcoup there. In case of partial kinematic coupling, velocities at the slab surface are low at the subduction channel end and increase downwards. z partial-coup is defined as the shallowest depth where velocities at slab surface are low but different from zero, while zcoup is shifted deeper at the level where overlying mantle velocities reach the subduction speed. C. The weak layer geometry imposed at simulation start is displayed in dark blue. The temperature field constantly maintained in the light blue dashed area simulates the lithosphere cooling from its formation at the ridge (top left corner) until a chosen lithosphere age, A lith , reached 400 km far away. A vertical resistance against flow is modelled along the open bottom box boundary, by mimicking a viscosity jump 10 times higher below [START_REF] Ribe | Three-dimensional modeling of plumelithosphere interaction[END_REF][START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF]: τ d zz is the deviatoric vertical stress, uz is the vertical velocity, and ν BB the viscosity at box bottom. 3), performed with the coarsest mesh LR (Table 1) and simulation S1-10b performed with the finest mesh HRlim. Depth (km) 3c1 and with a constant strain rate of 10 -14 s -1 . "B&K80": [START_REF] Brace | Limits on lithospheric stress imposed by laboratory experiments[END_REF], "Kal95": [START_REF] Kohlstedt | Strength of the lithosphere: Constraints imposed by laboratory experiments[END_REF], "R&M87": [START_REF] Ranalli | Rheological stratification of the lithosphere[END_REF], "K83": [START_REF] Kirby | Rheology of the lithosphere[END_REF], "Mal98": [START_REF] Mackwell | High-temperature deformation of dry diabase with applications to tectonics on Venus[END_REF], "W&C90": [START_REF] Wilks | Rheology of some continental lower crustal rocks[END_REF], "G&T95": [START_REF] Gleason | A flow law for dislocation creep of quartz aggregates determined with the molten salt cell[END_REF], "R95": [START_REF] Ranalli | Rheology of the Earth[END_REF], "Kal86": [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF]. (b) Down-dip extent of the interplate channel, zcoup, as a function of ∆Ea, for At set to 100 Myr and 60 Myr. The ∆Ea range for which partial kinematic couling is modelled is depicted in grey, and the range over which zcoup and z partial-coup vary through time by dashed and solid lines, respectively. In hachured areas, zcoup is mesh-dependent and located at 250 km depth, yielding an excessively cold fore-arc structure. 3). In the inset, the velocity color scale is saturated to 0.65 cm/yr. Right (b2): ∆Ea = 155kJ/mol (simulation S1-15). (b3): same simulation but including water transfer computation. The pink line is the 980 o C isotherm, i.e., the lowest solidus temperature for a water-saturated peridotite. "Eclg": slab dehydration triggered by the gabbroic crust eclogitization. "Chl out": water expelling by chlorite destabilisation within the initially serpentinised mantle layer. "Stable phase A": water retention in the subducting mantle as a result of phase A formation when serpentine breaks down. In the bottom inset, the logarithmic water content color scale highlights very low water amounts in the mantle wedge. 3). The boundaries between the different subduction regimes are arbitrarily extended by dashed lines (unexplored domains). The red coloured shadding stands for the fore-arc heating by corner flow decreasing as ∆Ea is increased.

σ b = C 0 + γ c ρ g z ε ref =10

Figure 1 :

 1 Figure1: Sketch of the subduction interplate and definitions of the main terms used to describe the mantle wedge structure (in brown). Green ascending arrows represent fluids expelled by the dehydrating slab. The red area depicts the mantle wedge part hot enough for melting in wet conditions. In the bottom panel, two kinds of kinematic coupling between the subducting slab and the overlying mantle are considered: In case of full coupling, velocities at the slab surface are close to the subduction speed as soon as the down-dip extent of the subduction channel is crossed, defining zcoup there. In case of partial kinematic coupling, velocities at the slab surface are low at the subduction channel end and increase downwards. z partial-coup is defined as the shallowest depth where velocities at slab surface are low but different from zero, while zcoup is shifted deeper at the level where overlying mantle velocities reach the subduction speed.

Figure 2 :

 2 Figure2: Kinematic and thermal boundary conditions and initial thermal state. The subducting plate velocity, v sub , is imposed on a 832 km wide and 13 km deep domain (green box), counted from the box left-hand side and from the box surface. One isotherm every 400 o C. The weak layer geometry imposed at simulation start is displayed in dark blue. The temperature field constantly maintained in the light blue dashed area simulates the lithosphere cooling from its formation at the ridge (top left corner) until a chosen lithosphere age, A lith , reached 400 km far away. A vertical resistance against flow is modelled along the open bottom box boundary, by mimicking a viscosity jump 10 times higher below[START_REF] Ribe | Three-dimensional modeling of plumelithosphere interaction[END_REF][START_REF] Arcay | Dynamics of interplate domain in subduction zones: influence of rheological parameters and subducting plate age[END_REF]: τ d zz is the deviatoric vertical stress, uz is the vertical velocity, and ν BB the viscosity at box bottom.

Figure 3 :

 3 Figure 3: Close-up on the subduction structure obtained in simulation S1-8 (Table3). (a) Thermal state (one white isotherm every 200 o C), composition (oceanic crust outlined by the green contour) and velocity field at different time steps. The trench is located at the abscissa x = 0 km. (b) Thermal structure and strain rate distribution at steady-state (26 Myr of convergence). One red outline every 8 × 10 -14 s -1 . (c): Close-up on the subduction interface. The interplate sampling line (orange line) joins points of maximum dissipation energy rates (see the text for details). Outlines of dissipation energy rate (purple lines) are depicted every 20 µW.m -3 . Panel(d): Profiles interpolated along the interplate sampling line of, from left to right, (1) second invariant of strain rate tensor, (2) temperature, and (3) viscosity.

Figure 4 :

 4 Figure 4: Influence of grid resolution and weak layer thickness on the down-dip extent of the decoupling interface: Zooms on the subduction structure at steady-state. Isotherms (black lines) every 200 o C. The rectangular Eulerian mesh is coloured as a function of crust content at computational nodes, by highlighting very high contents in crust material (>99.0 %). Within the white line the content in weak oceanic crust at Eulerian nodes is greater than 49.0 %. Panel (a): Simulation S1-10a (Table3), performed with the coarsest mesh LR (Table1) and simulation S1-10b performed with the finest mesh HRlim. Panel (b) Simulation S1-10c, computed with the finest mesh HR. Panel (c): Simulation S1-10d, performed with mesh MHR of intermediate resolution.

  Figure 4: Influence of grid resolution and weak layer thickness on the down-dip extent of the decoupling interface: Zooms on the subduction structure at steady-state. Isotherms (black lines) every 200 o C. The rectangular Eulerian mesh is coloured as a function of crust content at computational nodes, by highlighting very high contents in crust material (>99.0 %). Within the white line the content in weak oceanic crust at Eulerian nodes is greater than 49.0 %. Panel (a): Simulation S1-10a (Table3), performed with the coarsest mesh LR (Table1) and simulation S1-10b performed with the finest mesh HRlim. Panel (b) Simulation S1-10c, computed with the finest mesh HR. Panel (c): Simulation S1-10d, performed with mesh MHR of intermediate resolution.

Figure 6 :

 6 Figure 6: Influence of the crustal strength reduction, ∆Ea, on subduction mode. One black isotherm every 200 o C. The crust content at Eulerian nodes is >49.0% within the white outline. The numerical mesh is colored as a function of crust content, by focusing on very high contents (>99 %). Panel (a): Subduction mode with At =20 Myr depending on ∆Ea if the oceanic crust is very buoyant (∆ρc = -380 kg.m -3 , a1: simulation S2-8a, and a2: simulation S2-15a). Panel (b): Subduction of an old lithosphere for high ∆Ea. Stream function (green outlines) is depicted every 0.2 cm 2 /yr counted from -1.7 cm 2 /yr. Left (b1): ∆Ea = 145 kJ/mol (simulation S1-14, Table3). In the inset, the velocity color scale is saturated to 0.65 cm/yr. Right (b2): ∆Ea = 155kJ/mol (simulation S1-15). (b3): same simulation but including water transfer computation. The pink line is the 980 o C isotherm, i.e., the lowest solidus temperature for a water-saturated peridotite. "Eclg": slab dehydration triggered by the gabbroic crust eclogitization. "Chl out": water expelling by chlorite destabilisation within the initially serpentinised mantle layer. "Stable phase A": water retention in the subducting mantle as a result of phase A formation when serpentine breaks down. In the bottom inset, the logarithmic water content color scale highlights very low water amounts in the mantle wedge.

Figure 7 :

 7 Figure 7: Subduction modes for a 20 Myr old incoming lithosphere as a function of the subducting crust buoyancy, ∆ρc, and crust ductile weakening, ∆Ea. Experiments are depicted by small solid black circles (Table 3). In insets (1) to (7) illustrating the different subduction regimes, isotherms are displayed every 200 o C, and the subducting crust by the grey-dashed area. (1) simulation S2-12a in Table 3; (2) simulation S2-15c; (3) simulation S2-22c; (4) simulation S2-10h; (5) simulation S2-20j; (6) simulation S2-2c; (7) simulation S2-6g.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 891011 Figure 8: Evolution of the velocity field and thermal structure for a buoyant and weak subducting crust, when crust plumes form (panel a, simulation S2-13g, Table3), or not (panel b, simulation S2-17b). One isotherm (white outline) every 200 o C. Within the thick black outline, crust content is higher than 49%. The color scale refers to velocity magnitude, arrows depict velocity vectors.

Figure 12 :

 12 Figure 12: Subduction mode as a function of the incoming plate age and the contrast in ductile strength between crust and mantle, for a buoyant subducting crust (∆ρc = -380 kJ/mol) and a 100 km thick overriding plate. Numerical experiments are depicted by solid black circles (Table3). The boundaries between the different subduction regimes are arbitrarily extended by dashed lines (unexplored domains). The red coloured shadding stands for the fore-arc heating by corner flow decreasing as ∆Ea is increased.

Table 1 :

 1 Characteristics of numerical grids. for A t = 60 Myr and A t = 100 Myr, ∆z min is imposed for depth shallower than 250 km when ∆E a > 120 kJ/mol.

	Grid	LR	HRlim	HR	MHR
	Resolution x × z	332×90	482×129	482×151	407×119
	∆x min (km)	2.846	1.267	1.267	1.72
	imposed between (km) 1012 < x <1397 1004< x <1382 1004 < x < 1382 996 < x < 1389
	∆x max (km)	9.51	9.51	10.4	10.6
	∆z min (km)	2.268	1.266	1.266	1.656
	imposed for (km )	z > -116.1	z > -102.0	z > -136	z > -142 a
	∆z max (km)	10.2	10.3	10.3	10.1
	Oceanic crust
	thickness (km)	6.94	6.94	6.94	8.33
	Total tracer number	1.2544×10 6	5.76×10 6	5.76×10 6	3.24×10 6
	a

Table 2 :

 2 Parameter names and values. Density offset between crust and mantle ∆ρ c ∆ρ c = ρ c -ρ m Coupling depth at steady-state (t > 8 Myr) if A t > 20 Myr; z coup at 10 Myr of convergence if A t = 20Myr. Run includes the computation of passive water transfers. Run includes water transfer computation and a mantle strength reduction depending on water content.

	Parameter name	Symbol	Value
	Box height	H 0	555 km
	Bottom temperature	T b	1888 K
	Surface temperature	T s	273 K
	Mantle density	ρ m	3300 kg.m -3
	Crust density	ρ c	varying
	Mantle radiogenic heat production	A	9.20 ×10 -8 W.m -3
	Adiabatic gradient Thermal diffusivity	∂T ∂z adiab κ	0.45 K.km -1 0.8 ×10 -6 m 2 .s -1
	Thermal expansion coefficient	α	3.5 ×10 -5 K -1
	Heat capacity	C p	0.971 ×10 3 J.(K.kg) -1
	Activation energy for the mantle Activation energy for the crust Activation energy offset crust/mantle Activation volume	E m a E c a ∆E a V a	465 kJ/mol varying ∆E a = E m a -E c a 1.7×10 -5 m 3 /mol
	Pre-exponential factor in non-Newtonian rheology	A 0	339428.7 MPa -3 .s -1
	Dissipation number	Di = αgH0 Cp	0.196
	Gravity acceleration	g	9.81 m.s -2
	Weak layer thickness	h c	6.9 km or 8.3 km
	Convergence rate	v c	6.5 cm/yr
	Subducting plate age at trench	A t	varying
	Cohesive strength at z = 0	τ 0	1 MPa
	Stress exponent in the viscous rheology	n	3
	Stress exponent in the brittle rheology	n p	30
	Yield stress increase with depth (mantle)	γ m	1.6
	Yield stress increase with depth (crust)	γ c	0.065

b Coupling depth dependent on numerical mesh. c d Flat and aborted subduction.

e

Acknowledgement

Comments from Satoru Honda, James Conder, an anonymous reviewer, and editor Marc Jellinek have lead to a significant improvement of the paper. Their careful reviews were greatly appreciated. Stéphane Arnal, Fabrice Grosbeau, and Josiane Tack are acknowledged for the maintenance and development of the cluster of computing nodes on which all numerical experiments were performed. I deeply thank Matthieu Lilla, Izia L. and Elio L., as well as Marc Regny-Demery for their support that was essential. Figure 9 is made using GMT (Wessel and Smith, 1991). This work has been funded by the CNRS-INSU

. Lines with stars depict the hottest geotherm modelled in the mantle wedge, compared with conditions of last equilibration with the mantle of near-primary arc magmas [START_REF] Tatsumi | Melting experiments on a high-magnesium andesite[END_REF][START_REF] Tatsumi | Generation of arc baslt magmas and thermal struction of the mantle wedge in subduction zones[END_REF][START_REF] Bartels | High pressure phase relations of primitive high-alumina basalts from Medicine Lake Volcano, Northern California[END_REF][START_REF] Draper | Anhydrous P-T phase relations of an Aleutian high-MgO basalt: An investigation of the role of olivine-liquid reaction in the generation of arc high-alumina basalts[END_REF][START_REF] Falloon | Melting of refractory mantle at 1[END_REF][START_REF] Elkins Tanton | Hot, shallow melting under the Cascades volcanic arc[END_REF][START_REF] Grove | Fractional crystalization and mantle-melting controls on calcalkaline differenciation trends[END_REF][START_REF] Hesse | Absakorites from the western Mexican vocalnic belt: constraints on mantle wedge conditions[END_REF][START_REF] Watt | Arc magma compositions controlled by linked thermal and chemical gradients above the subducting slab[END_REF], red solid squares, red-dotted box).