Interactive comment on “ Dynamics of interplate domain in subduction zones : influence of rheological parameters and subducting plate age ” by D

I thank Referee 2 for her very constructive comments. I agree that results were initially badly illustrated by the choice of figures. I follow the clever suggestions made by Referee 2 and build a plot of the modelled zBDT obtained as a function of zdec in simulations of old subducting lithosphere subduction (new figure, with label 5). It is then much easier to illustrate the respective influence of each rheological parameters on zBDT and zdec. I recognise that former Figure 9 was awkward. I replace it by the same type of plot as mentioned just above, where zBDT is presented versus zdec for both young and old subducting plates. This highlights the influence of lithosphere


Introduction
The subduction interplate domain (considered either as a plane or a channel, depending on the setting) is an interface of seismogenic coupling at the time scale of one seismic cycle, and also of kinematic decoupling on long-term geological time scales.The properties of this very particular interface are likely to affect not only the seismogenic potential of the subduction area but also the overall subduction process, as it influences its viability.However, the different mechanisms governing the subduction interplate dynamics remain poorly known.For instance, a great variability of the down-dip limit of the seismogenic zone is observed, encompassed between 30 and 70 km (e.g.Pacheco et al., 1993;Heuret et al., 2011).This depth might be (at least partly) controlled by the brittle-ductile transition occurring along The light-dotted domain depicts the motionless upper lithosphere.The subduction interplate plane, here envisioned as a tangential kinematic discontinuity, is coloured in black.The interplate weak layer, located at the subducting lithosphere surface as an oceanic subducting crust layer, is depicted in grey.The "viscous blanket" refers to the thermal boundary layer formed by asthenospheric cooling at the subducting lithosphere surface (Kincaid and Sacks, 1997).(b) Definition of the brittle-ductile transition depth, z BDT .The brittle-ductile boundary (red line) connects rock elements where the pseudo-brittle strength, ν b , equals the non-Newtonian viscosity, ν v .The medium is modelled as brittle above and ductile below, as sketch in the stress-depth diagram along the interplate plane on the right.The shallowing effect on z BDT of an energy activation decrease, keeping constant the reference viscosity, is sketched in blue.
the subduction channel, and could thus depend on many variables such as temperature, pressure, compositional variations, strain rate, etc.This implies a self-consistent equilibrium state of the subduction interplate, whose characteristics would depend on the subduction setting.Numerical simulations of subduction dynamics appear as one of the more powerful tools to try to unravel the physics of the interplate dynamics.

Modelling the subduction interplate in simulations of convergence
Numerical modelling of subduction shows that the technique used to simulate the kinematic decoupling between the two converging plates has a huge influence on the produced features interesting the modeller, especially at the convective mantle wedge tip where the contact between the two plates stops (Fig. 1a).There occur very high gradients in temperature, strain rate, and strength that eventually govern the most characteristic patterns of the mantle wedge dynamics detected at the surface: heat flow increase in a domain of overall cooling, partial melting, and high flux of expelled fluids.Also, numerical models demonstrate that the low strength fault zone decoupling mechanically the two plates has to be assigned to mimic a realistic convergence zone; otherwise, a complete viscous mechanical coupling between plates takes place and the lower part of the fore-arc mantle is ablated, yielding an extreme heating at shallow depth (e.g.Eberle et al., 2002).Several methods have been explored to decouple kinematically the two converging plates.A first set is based on kinematic assumptions, such as imposing: free slip along the boundary (Furukawa, 1993); rigid and motionless fore-arc lithosphere (Peacock and Hyndman, 1999;van Keken et al., 2002), or a progressive kinematic coupling between the upper lithosphere sublayer and the subducting slab (Kneller et al., 2005(Kneller et al., , 2007;;Syracuse et al., 2010).Another approach aims at simulating low strength/low shear along the interplate boundary, by either assigning low viscosities to the interface nodes (Billen and Gurnis, 2001;Kelemen et al., 2003;Wada et al., 2008;Wada and Wang, 2009), or limiting shear stress (Zhong and Gurnis, 1995;van Hunen et al., 2002), or impeding fore-arc deformation if predicted to occur in the brittle domain, itself being delimited by a predetermined temperature (Conder, 2005;Syracuse et al., 2010).Thermo-kinematic models with prescribed interplate mechanics are useful to test specific assumptions suggested by observations, such as partial melting domain extent and/or geometry of the cold fore-arc nose, because subduction geometry can easily be adjusted to fit the observed one.However, one part of the involved physics regulating the interplate equilibrium cannot be resolved and demands a dynamic modelling in which temperature, flow, and stress evolve freely and consistently as a function of their own interactions.This is the main purpose motivating this paper.
In the present study, numerical models are performed to study the equilibrium state of the subduction interplate when the diving lithosphere interacts with both the overriding plate and the surrounding convective mantle, after a 650-900 km length of subduction, i.e. when the subduction transient state (more or less) ends.The decoupling interface geometry is not fixed and its properties are not assigned, as both evolve as a function of advection of weak crustal material within the interplate layer.Rock strength here depends on depth, temperature and stress, for both the mantle and the weak crust filling the interplate domain.The thermo-mechanical model combines a non-Newtonian viscous rheology and a pseudo-brittle rheology.By combining these two mechanical behaviours, one is then able to study how the bounds of the brittle realm along the subduction plane, on the one hand, and the down-dip extent of kinematic decoupling between the two converging lithospheres, on the other hand, stabilise through time and possibly interact, as a function of (1) rheological parameters and (2) subduction parameters (e.g.convergence rate, subducting lithosphere and upper plate structures, asthenosphere flows).Regarding item 2, this paper focuses on the influence of the subducting lithosphere age.

Depth of interplate kinematic decoupling vs. depth of brittle-ductile transition
Kinematically speaking, the tangential displacement between the upper lithosphere and the subducting slab is decoupled on both sides of the interplate plane.Below the interplate decoupling depth, mantle rocks overlying the subducting slab are passively dragged down by the latter.The transition depth between decoupled motions above and coupled displacements below is labelled the "interplate decoupling depth" (Furukawa, 1993).Advection of warm asthenospheric rocks occurs into the wedge to replace the mantle dragged down along the slab by viscous coupling across the slab top (labelled corner flow, Fig. 1a).This rising return flow is mainly passive.As a result, a large temperature jump occurs across the slab top in the vicinity of the interplate decoupling depth, resulting in a drastic mantle viscosity decrease if the rheology is non-Newtonian and temperaturedependent (Andrews and Sleep, 1974;Honda, 1985), and entails a corner flow focussing at the decoupling interface base.Moreover, focused high strain rates, confined in the decoupling interface until its down-dip extent, jump away from the slab surface and reach the asthenospheric wedge over a relatively narrow interval where thermal gradients are very high (Fig. 1a).Therefore, the interplate decoupling depth results from a thermomechanical equilibrium, probably depending on the asthenosphere/interplate material strength contrast, and also on subduction parameters, such as convergence rate, that govern the interplate strain rates and flow velocities in the mantle wedge.
From the surface to a given depth along the interplate plane, stress along the subduction plane increases with depth together with the brittle yield stress.Meanwhile, temperature increases and reduces the ductile strength.As a consequence, the brittle stress increase finally crosses the ductile stress decrease at a depth where interplate stress is maximum, defining the brittle-ductile transition (BDT, Fig. 1b).The depth of BDT, z BDT , cannot exceed the interplate decoupling depth, z dec , but some authors have assumed that the kinematic decoupling was occurring at the BDT (Furukawa, 1993;Conder, 2005;Arcay et al., 2007bArcay et al., ,a, 2008)).I test this hypothesis in this paper and show that conditions required to simulate z BDT = z dec may be much more restricted than initially expected.
Note that the BDT may result from a complex phenomenon, possibly involving metamorphic reactions and fluid migration, as suggested by non-volcanic tremors (e.g.Obara, 2002;Audet et al., 2009) and/or specific deformation mechanisms (e.g.Branlund et al., 2001), which will not be tested here.The interplate domain is here simply modelled by a layer compositionally different from the surrounding mantle, with specific rheological parameters.In nature this interface is probably made of pounded material mixing subducted sediments and slices of over-riding crust torn during underthrusting (e.g.Lallemand et al., 1992;Lallemand, 1995;Collot et al., 2011), therefore much weaker than the subducting oceanic crust.Hence, its rheological properties are assumed to be close to those of a continental crust to mimic the behaviour of a real subduction channel.From a technical point of view, it is nevertheless easier to assume that the layer localising deformation has the geometry of the oceanic subducting crust.Its density must however be adjusted as if it were oceanic crust to model correctly the slab pull and a realistic force balance.
The paper starts with the description of the modelling setup.Next, the dynamics of the subduction interplate is simulated for two end-member ages of the subducting lithosphere, 100 Myr and ∼ 20 Myr, representing the interval encountered on Earth in the vast majority of subduction zones (Heuret and Lallemand, 2005) (excluding three subduction zones: Cascadia, Mexico, and the Chile triple junction where the subducting plate is younger than 15 Myr).I first study how z BDT and z dec equilibrate for a subducting plate ∼ 100 km thick, by varying brittle parameters and also the non-Newtonian strength.I then lean on the derived conclusions to explain why the subduction of a young lithosphere may be sustainable or not, depending on the modelled rheology.The influence of the convergence rate on z dec has been already extensively studied elsewhere (Arcay et al., 2007b(Arcay et al., ,a, 2008) ) and is not investigated here.

Model setup
A thermochemical code of convection (Christensen and Yuen, 1984;Christensen, 1992) is used to model subduction.It solves the momentum, energy, and mass conservation equations.Rocks are assumed to be incompressible, except for the thermal buoyancy term in the momentum equation, and for the adiabatic heating term (Table 1) in the energy equation (extended Boussinesq approximation).Shear heating (i.e.viscous and frictional dissipations) and uniform heat production are also included in the heat conservation equation.Indeed, shear heating has been shown to help significantly strain localisation and weak strength inside the subduction interface by sustaining high temperature (e.g.Doin and Henry, 2001;Faccenda et al., 2008) box is 2220 km wide and 555 km high (Fig. 2).Composition (either mantle or weak crust) is tracked by two types of tracer that have different densities and rheological properties.Buoyancy depends on temperature, through the thermal expansion coefficient, and composition (crust/mantle).Compositional tracers are advected with the velocity field (van Keken et al., 1997).

Mechanical boundary conditions and subduction modelling
Subduction is simulated by applying a constant convergence rate of v sub = 6.5 cm yr −1 on top of the incoming lithosphere, on a 832-km-wide and 16-km-deep segment (Fig. 2).The diving plate then evolves freely within the trench area.The upper lithosphere is 100 km thick and is here assumed to be simply made of mantle rocks.The incoming plate is covered by a 7-km-thick layer of "crust" material much weaker than the underlying mantle.At simulation start, an initial 30 • dipping interplate layer made of weak crust material is imposed from the surface to 55 km depth, at the middle of the box.The trench is hence initially located 1110 km away from the lefthand side of the box.Strain localisation along the interplate boundary during convergence basically relies on the strength contrast between the weak layer plane and the mantle, composing the upper lithosphere and also the remaining part of the subducting plate.Deformation localisation along the convergence boundary is then a function not only of the specific mechanical properties of the modelled mantle and crust, but also on the interplate geotherm or, more precisely, on the difference between interplate and upper plate geotherms.If the thermo-mechanical conditions are such that the weak crustal material is able to localise deformation, the subducting mantle lithosphere bends, and the weak material flows at the subducting plate surface to continuously fill the initial interplate channel.Subduction in this case is successfully initiated and is sustained by the constant convergence rate.
The lower box boundary is open to prevent unrealistic slab deformation that would occur if the slab encountered the box base.However, a vertical resistance against flow is modelled in some simulations to help the convective mantle, if not resistant enough, to compensate the subducting slab weight.If k is the wavenumber of a harmonic vertical flow field, the resistance to vertical flows applied at box bottom, R b , writes as R b = ν × k, where ν is the viscosity of the virtual material that should underlie the open lower boundary.The boundary condition along the box bottom is σ zz −2R b v z = 0, where σ zz is the non-hydrostatic vertical stress and v z is the vertical velocity (Ribe and Christensen, 1994).By setting the reference strain rate in the model, εref , to 10 −14 s −1 , and the subduction velocity, v sub , to 6.5 cm yr −1 , the vertical scale length of deformation is L = v sub /ε ref , which defines here the main wavelength of deformation, k −1 .ν is set to either the normal viscosity at box bottom, ν BB : ν * BB = 1 in ν = ν BB × ν * BB (no viscosity jump across the lower boundary, simulations S10 and S12 for instance in Table 2) or 10 times the viscosity at 555 km depth (ν * BB = 10, e.g.simulations S12 and S13, Table 2).Other mechanical boundary conditions are presented in Fig. 2.

Thermal boundary conditions -modelling the subduction of a constant age-at-trench lithosphere
The whole convective box is heated by a uniform radiogenic source (Table 1).Along the surface, the temperature is set to 0 • C, whereas all other boundaries are insulating.At subduction initiation, the upper lithosphere thermal structure is that of an old ( 100 Myr) and cold lithosphere, at equilibrium with the underlying convective mantle, thence stable.The same thermal state is applied for the incoming plate in models of 100-Myr-old lithosphere subduction (section 3).This thermal structure ensures constant equilibrium between cooling from above and heating from below by asthenosphere convective flows, which prevent any plate thermal thickening during subduction.The simulation box thermal structure is the result of a preliminary run in the same condi-tions as described above, but without convergence velocity.
The lithosphere has finally a homogeneous 100 km thickness (Fig. 2a), with basal small perturbations resulting from smallscale convection.A pseudo-ridge is imposed at the incoming lithosphere extremity, and simulates the plate conductive cooling from 0 Ma to a chosen lithosphere age, A lith , on a 400 km width (Fig. 2).The structure of this pseudo-ridge is constantly sustained as a boundary condition, and is used to re-generate the incoming lithosphere while it is consumed by subduction.An overlying layer of 7-km, weak, crustal material is also constantly maintained on the surface of the newly formed lithosphere.As imposed by the assigned boundary conditions, a segment of lithosphere of age A lith is located 710 km away from the trench and undergoes an ageing of 11 My during its route to the subduction trench at a speed of 6.5 cm yr −1 .As a consequence, the value of A lith is set to 90 Myr to account for the newly formed lithosphere cooling and thickening and to finally model a subducting lithosphere of constant age.
The thermal thickness of a 20-Myr-old subducting lithosphere, defined by the 1200 • C isotherm depth, should be close to 52 km according to the half-space cooling model (Turcotte and Schubert, 1982).However, this thickness results in a predicted surface heat flux of 69 mW m −2 for the thermal conductivity I assume (Table 1), while surface heat flux estimates of oceanic basin floor indicate a value of rather 112 mW m −2 for a 20-Myr-old oceanic lithosphere (Doin and Fleitout, 1996a).This high heat flux would imply a quite thin lithosphere only 31 km thick.A compromise temperature gradient of 34.6 • C km −1 is imposed from the   Solid Earth, 3, 467-488, 2012 www.solid-earth.net/3/467/2012/surface to the lithosphere base (39 km depth) to mimic a 20-Myr-old plate.A lith is then adjusted to 10 Myr to maintain a roughly constant lithospheric age of 20 Myr at trench.Note that the thermal structure of the interplate area is a bit cooled at subduction onset for hot incoming plates to help strain localisation.The isotherms of the incoming lithosphere are curved to be parallel to the subduction plane, which enhances the strength of the first subducting lithospheric segment and favours deformation localisation within the weak plane (Fig. 2b).

Rheological model
The rheological model combines a pseudo-brittle rheology, with a yield stress increasing with depth, to a non-Newtonian creep rheology.An effective viscosity, ν eff , is defined through the relationship: τ = ν eff ε, where τ and ε are the second invariants of the stress and strain rate tensors, respectively.The effective viscosity is given by the harmonic average between a brittle-plastic term, ν b , and the non-Newtonian viscosity strength, , assuming that the total deformation is the sum of brittle and ductile strains.
The pseudo-brittle rheology is modelled through a yield stress, τ y , increasing with depth, z: τ y = τ 0 + γ (C)ρgz, where τ is the cohesive strength at the surface and γ is a coefficient depending here on composition.This coefficient is related to the friction coefficient, f s , through the relationship (Turcotte and Schubert, 1982) , where ρ w is the water density.The effective plastic viscosity, ν b , is given by where εref is a reference strain rate and n p is a large stress coefficient (Table 1).In the plastic domain, very large strain rates are simulated as soon as stress exceeds the yield stress.
The non-Newtonian rock viscosity, ν v , writes as where T is temperature in kelvin, A 0 the pre-exponential constant, E a the activation energy, depending on composition, V a the activation volume, n an exponent different from 1, and R is the gas constant (Table 1).The choice of a non-Newtonian rheology in the creep deformation model favours the development of a self-sustaining localisation of deformation within area of high strain-rate and low strength and facilitates subduction initiation (Billen and Hirth, 2005).
In two simulations (S13f14 and S13f14b, Table 2), a local viscosity decrease is modelled in the mantle wedge area by simulating (1) slab dehydration related to destabilisation of oceanic crust hydrous minerals and deserpentinisation of the underlying subducting cold mantle, (2) vertical migration of expelled fluids and subsequent water-saturation of the overlying asthenospheric rocks in the hydrated mantle wedge area, and (3) hydrous strength weakening for hydrated mantle rocks.The model of slab dehydration prediction-mantle wedge hydration in the absence of slab melting based on accurate P − T phase diagrams was extensively described in Arcay et al. (2005Arcay et al. ( , 2006)).I do not want to go into the details of the modelling, since the purpose in this paper is to use it as a way to simulate a viscosity decrease only in the vicinity of the mantle wedge tip where the subduction interface ends.In a nutshell, the dehydration-hydration geometry is basically a function of the diving lithosphere thermal state.For a 6.5 cm yr −1 convergence rate and a 100-Myr-old plate, the hydrated area width in the mantle wedge equals at maximum 133 km, 6.5 Myr after subduction initiation, and decreases to the steady value of 84 km at 15 Myr, as a consequence of the subduction dip angle increase while slab pull develops.I choose to reduce for water-saturated rocks the pre-exponential constant in Eq. 2 to A 0 /14, which results in a hydrous strength reduction of 14 3/2 ∼ 52 (Table 2) if the effective non-Newtonian viscosity has to be expressed rather as a function of the energy dissipation rate, (σ ε), than strain rate, ε (Christensen, 1984;Dumoulin et al., 1999).

Choice of rheological parameter sets
The goal is to investigate the interplay between brittle and ductile rheologies along the subduction interface.As brittle deformation is mainly controlled by the friction coefficient, the only tested parameter is the weak layer frictional coefficient, γ c , keeping constant the cohesive strength at the surface, C 0 , the exponent n p in Eq. 1 and the reference strain rate, εref .To simplify, the mantle friction coefficient, γ m , will also remain set to 1.6 (Table 1).
In this paper, creep behaviour is a function of temperature, pressure, strain rate, and composition (mantle/weak layer).The temperature-and pressure-dependence in viscosity are respectively controlled by the activation energy and the activation volume (Eq.2).The effects of V a , weak crust and mantle activation energies, E c a and E m a , have to be tested.Activation volume and activation energy represent only strength gradients in the logarithm of viscosity (associated with thermal and pressure gradients, respectively), and a reference viscosity at a given depth and temperature has to be assigned.It has been shown that the effective viscosity at the lithosphere base, corresponding to a minimum if T -and P -dependences are both modelled, is crucial in triggering of small-scale convective processes (e.g.Davaille and Jaupart, 1993;Dumoulin et al., 2001;Morency et al., 2002;Dumoulin et al., 2005), as well as in subduction initiation, as it favours, if low enough, the corner flow activation necessary to stop the kinematic decoupling and mechanical coupling between the two converging plates at high depths (e.g.Doin www.solid-earth.net/3/467/2012/Solid Earth, 3, 467-488, 2012 and Henry, 2001; Kukačka andMatyska, 2004, 2008).Therefore, I choose to scale ductile strengths with different activation energies and volumes by adjusting the pre-exponential constant, A 0 , in Eq. ( 2), to keep the asthenosphere viscosity at the lithosphere base, ν asth (z =100 km depth, T = 1350 • C, εref = 10 −14 s −1 ), equal to 2.724×10 19 Pa•s.Nevertheless, the influence of ν asth is also tested by modelling in a few experiments a hydrous strength weakening associated with mantle wedge metasomatism (see Sect. 2.3.1).Hence, four different ductile rheologies are tested (labelled C6, C10, C12, and C13, Table 2), in order to investigate the respective influence of weak crust activation energy, E c a , mantle activation energy, E m a , the difference between the latter two, E a = E m a − Ea c , and the activation volume of both compositions, V a .Rheology C6 simulates a weak crust whose strength is close to the undried Adirondack granulite studied by Wilks and Carter (1990) and a mantle strength similar to the one of wet synthetic olivine (Karato et al., 1986).Rheologies C6 and C12 are identical except for a moderate difference in activation volume, and both simulate the ductile behaviours of, respectively, the dry mafic granulite (Wilks and Carter, 1990) for the weak layer, and an intermediate strength between wet dunite at low temperature (Chopra and Paterson, 1981) and wet olivine (Hirth and Kohlstedt, 1996) at high temperature.Rheology C13 models a weak crust close to the dry diabase studied by Kirby (1983) and a mantle with a strength similar to wet Aheim dunite at high temperature (Wilks and Carter, 1990) and encompassed at low temperature between wet olivine in Kirby (1983) and wet synthetic olivine in Karato et al. (1986).
Finally, since the depth of brittle-ductile transition may depend on interactions between ductile creep and brittle deformation, two strategies are followed to study this possible interplay.On the one hand, the ductile thermo-mechanical parameters are kept constant, whereas I vary the crust frictional coefficient (rheologies C10 vs. C10LG; C13 vs. C13HG; C13 vs. C14b, Table 2).On the other hand, the friction parameter is kept constant in rheologies C6 and C12 (and in C13 and C13fnu14), whereas the combination of ductile rheology parameters is changed.

Numerical resolution
The conservation equations are solved with a spline finite element method on a non-deforming grid (Eulerian approach, Christensen and Yuen, 1984;Christensen, 1992).The simulation box, 2220 km wide and 555 km high, is discretised into 332 × 90 nodes.The grid is refined to improve resolution in areas of high thermal and deformation gradients.Close to the subduction plane and in the wedge tip area, the horizontal and depth grid spacings are 2.8 and 2.3 km, respectively.Outside the mantle wedge, they are equal to 9.5 and 10.2 km, respectively.The tracer density is uniform over the simulation box (1 per km 2 ), with a minimum of seven tracers in the smallest meshes.The numerical discretisation used here as well as the tracer density was validated in a former study (Arcay et al., 2005).

Interplate dynamics for an old incoming lithosphere
The two following sections describe how the interplate structure evolves as a function of rheological parameters when the subducting lithosphere is old and ∼ 100 km thick, i.e. when the strength contrast between the weak interplate material and the surrounding mantle is maximum.The weak "crust" density is set to 2920 kg m −3 .The time window is encompassed between 10 and 14 Myr, i.e. close to the end of the subduction transient state, to capture both the main characteristics of the subduction steady-state and of their evolution.

Estimate of z BDT and z dec
Let us first define the brittle-ductile transition depth, z BDT , and the interplate decoupling depth, z dec , in these simulations.Figure 3a illustrates the subduction zone in simulation S13, 16 Myr after simulation start.As the subduction interplate domain concentrates deformation, it corresponds to a maximum in energy dissipation rate, compared to neighbouring areas (Fig. 3b).The coordinates of the subduction plane are thus extracted by looking for maximum dissipation rates in a window encompassing the converging boundary, and define the interplate sampling line depicted in Fig. 3a  and b.Four mechanical fields are then interpolated along this line: dissipation rates, deviatoric stress, strain rate, and temperature (Fig. 3c).Since brittle strength increases linearly from the surface to a maximum at the brittle-ductile transition depth, z BDT is defined by the depth of maximum stress (Fig. 3c2).I verify that this depth is never deeper than the brittle-ductile boundary defined by the relationship ν v (x, z) = ν b (x, z) (red line in Fig. 3b).Beyond the z BDT depth, deformation is mainly ductile and stress decreases while temperature keeps on rising.At the down-dip limit of the subduction plane, the slab surface is in contact with the asthenosphere and there occurs a strong increase in thermal gradient.Asthenospheric rocks become there mechanically coupled to the subducting plate.z dec is hence defined by the depth where the strongest temperature increase is recorded along the interplate sampling line, which should correspond to the location where isotherms at the mantle wedge tip become sub-vertical (Fig. 1a and 3b).As deformation jumps from the subduction interface to the asthenospheric tip at the depth z = z dec , we verify that strain rate, dissipation rate and stress decrease to zero for z > z dec (Fig. 3c).

Rheological parameters controlling the brittle-ductile transition depth
Figure 4a summarises the modelled z BDT for all tested rheologies obtained between 10 and 14 Myr after subduction onset.The brittle-ductile transition is usually not stabilised yet   4) temperature.Note that, in this particular case, the brittle-ductile transition depth, z BDT , is so close to the interplate decoupling depth, z dec , that they are assumed to be equal.
in most cases, and either deepens or shallows through time, but the order from the shallowest to the deepest BDT is generally constant.I focus on z BDT obtained 12 Myr after subduction initiation to discuss the effects of rheological parameters.
First, as one may expect, for constant asthenosphere viscosity at the down-dip limit of the subduction plane (i.e.aside rheology C13f14), the BDT depth is mostly dependent on the crustal friction coefficient, which defines the slope of yield strength increase (Fig. 1b).z BDT is minimum for the highest friction coefficient (z BDT =47.9-50.4km, γ c = 0.069, rheologies C10 and C14b) and progressively deepens as γ c is reduced in rheology C13HG (z BDT =51.8 ± 7 km, γ c = 0.061), C12 (z BDT = 58 km, γ c = 0.045), C6 (z BDT = 64 km, γ c = 0.045), C13 (z BDT = 77 km, γ c = 0.034), and C10LG (γ c = 0.007, z BDT ∼84 km but the measurement is quite inaccurate for this extremely low friction coefficient).When friction coefficient decreases, the brittle stress envelope intersects the ductile stress curve at deeper levels (Fig. 1b and Fig. 5, compare modelled brittle strengths in simulations S13 and S14b), assuming that the interplate ductile strength is roughly constant at a given depth from a simulation to the other.The only exception to this rule is obtained in simulation S14b, where z BDT is always deeper than in simulation S10 (at minimum by 2.6 km at 12 Myr, up to 11 km at 10 Myr), despite the same friction coefficient (γ c =0.069).The influence of the weak material activation energy appears here: if ductile strengths are scaled to keep the asthenosphere viscosity constant, a lower activation energy (E c a (C10) = 335 kJ mol −1 , while E c a (C14b) = 360 kJ mol −1 ) reduces the interplate viscosity decrease for a given geotherm, and yields a z BDT shallowing (Fig. 1b).Note that determining z BDT in simulation S14b is not straightforward (as in simulations S10 and S13HG), because stress does not strongly decrease after the maximum stress and slowly lowers with depth, until the abrupt reduction once z dec is exceeded (Fig. 5a).The depth  2).The question mark underlines the assumed estimate for z BDT at 12 Myr for rheology C10LG, with an extremely low friction coefficient for the weak layer, which makes the measurement impossible at that time.Minimum and maximum values can nevertheless be measured, as illustrated by the uncertainty bar.
The assumed z BDT value at 12 Myr is set to mid-depth in the uncertainty interval, as it also corresponds to the interplate decoupling depth at this time (b).It is indeed likely that z BDT = z dec in this case.(b) Interplate decoupling depth evolution, for an old subducting lithosphere, simulated for the same models as in (a).
interval between z BDT and z dec in which stress gradually reduces in simulations S10, S13HG, and S14b is interpreted as the crust ductile domain (Fig. 5c).
Apart from γ c , the BDT depth is strongly influenced by the asthenospheric viscosity.Simulations S13 and S13f14, differing by the mantle viscosity at the asthenospheric wedge tip, clearly show a z BDT shallowing (by ∼ 11 km) associated with the imposed mantle viscosity reduction.As depicted in Fig. 4b, the BDT depth equals the interplate decoupling depth at all times in simulations S13 and S13f14.The interpretation in that specific situation is that the BDT is controlled by the kinematic decoupling depth, z dec (see Sect. 3.4).A softer asthenosphere is more easily entrained by the slab surface, which shallows the kinematic coupling between wedge mantle and subducting plate and results in a z dec decrease.As a result, the enhanced heating at the subduction plate down-dip extent softens the weak layer ductile strength (Fig. 5, compare simulations S13 and S3f14), and finally brings about the transition from brittle to ductile deformation closer to the trench (Arcay et al., 2007b(Arcay et al., , 2008)).
At last, simulations show that an activation volume increase shallows z BDT (compare simulations S6, V a = 1.5 × 10 −5 m 3 mol −1 , and S12, V a = 1.7 × 10 −5 m 3 mol −1 ).Note that we obtain z BDT = z dec as soon as 10 Myr, but only from 12 Myr in simulation S6 (Fig. 4a and b).Therefore, the influence of V a on z BDT may be the direct consequence of the interplate decoupling depth dependence in V a (see Sect. 3.3).Fig. 6 summarises the influence of rheological parameters on z BDT and z dec for a 100-Myr-old subducting lithosphere.

Rheological parameters controlling the interplate decoupling depth
The simulated interplate decoupling depth is always encompassed in a depth interval narrower than obtained for the BDT (between 64 and 92 km depth, whereas the z BDT interval is 40-92 km, Fig. 4b).

Parameter with weak (to very weak) influence on z dec : crust friction coefficient
Considering models with constant asthenospheric viscosity, simulations C13, C13HG, and C14b model very close z dec values (∼87 km depth), despite very different friction coefficients (Table 2).Similarly, simulations C10 and C10LG with identical ductile rheological parameters but distinct friction coefficients converge towards the same z dec depth for subduction durations longer than 12 Myr, suggesting that the friction coefficient influence vanishes.Likewise, simulations S10 and S14b with the same crust friction coefficient but different ductile parameters show very different kinematic decoupling depths.I conclude that the interplate decoupling depth is thus independent from γ c .

Parameters governing z dec : asthenospheric strength, activation volume, and activation energy discrepancy between weak crust and mantle
By definition, z dec is the maximum depth of strain localisation along the interplate: for z > z dec , maximum strain rates jump to the viscous blanket surface and the upper plate base, but are not located within the subducting crust anymore (Fig. 1a).This is further highlighted by the interplate sampling line, materialising maximum dissipation rates along the interplate, that jumps away from the weak layer precisely at the depth z = z dec (Fig 3b).The fundamental role of the hot corner flow in subduction sustainability is to stop the kinematic decoupling between the converging structures (Doin and Henry, 2001): both the compressive setting and the thermal adherence between the cold subducting plate and the cooled fore-arc maintain the mechanical coupling between the converging plates until high depth and act to prevent a dipping and self-sustained subduction.This is the reason why the viscosity at the lithosphere-asthenosphere boundary has a significant effect on z dec .The hydrous viscosity reduction modelled in simulation S13f14 promotes the asthenosphere drag by the subducting slab, which in turn amplifies the heat flux advected by the return flow towards the interplate downdip extent.Similarly, as mentioned above, a slight increase in activation volume results in a moderate z dec shallowing (by ∼6 km) (Fig. 4b).As the viscosity at the upper plate base is identical in rheology C6 and C12, the V a increase emphasises the viscosity contrast with respect to shallower and deeper strengths, which enhances the corner flow focussing at the mantle wedge tip (Kukačka and Matyska, 2008), and shallows a bit the maximum depth of plate kinematic decoupling.
The last parameter affecting the interplate decoupling depth is neither the activation energy of the weak crust material, E c a , nor the mantle one, E m a , alone, but is the dif- Increase in crust activation energy Fig. 6. brittle-ductile transition depth versus interplate decoupling depth modelled after 12 Myr of convergence, for an old subducting lithosphere (simulations S6, S10, S10LG, S12, S13, S13HG, S13f14, and S14 in Table 2).Both depths are equal along the thick black line.The influence of rheological parameters on z BDT and z dec is highlighted in red and green, respectively.
ference between the latter two, E a , if mantle strength, ν asth , at the mantle wedge tip, is unchanged (simulations S6, S12, S13, and S10).z dec is minimum for the lowest difference in activation energy between mantle and weak crust ( E a = 65 kJ mol −1 , z dec =58 and 64 km, rheologies C12 and C6), deepens to ∼ 87 km when E a increases (105 kJ mol −1 , for rheologies C13, C13HG and C14b), until its deepest value modelled for the highest E a in rheology C10 (155 kJ mol −1 , z dec =92 km).To understand this result, note that the kinematic transition from decoupling to coupling between the mantle wedge and the subducting plate is controlled by the strength contrast (1) between the interplate domain and the upper lithosphere, and (2) between the subducting slab surface and the mantle wedge.Let us consider the interplate segment encompassed between z BDT and z dec in the ductile realm.The weak crustal material localises deformation along the subduction interface, because it is weaker than the surrounding cooled mantle.However, the interplate geotherm is much colder than a classical vertical geotherm across a 100-km-thick oceanic lithosphere (see Fig. 3c4).At the depth z = z dec , the cold subducted crustal layer is suddenly strongly heated but never reaches the asthenospheric temperature (the weak layer temperature being always lower than 800 • C, Fig. 3b).As a consequence, at the depth of the lithosphere-asthenosphere boundary, the weakest layer is the asthenosphere and not the subducted crust anymore, implying there is a jump of strain localisation towards the hot mantle (Fig. 1b).Of course, the contrast between mantle and crust ductile strengths is given by E a , meaning that a high E a promotes strain localisation within the weak subducted crust until deep levels, which deepens z dec .Note that the corner flow is active if the asthenosphere is successfully dragged downward by the slab, which occurs only if the dragging layer at the slab-asthenosphere contact is more viscous than the mantle wedge medium (as clearly evidenced by Kukačka and Matyska, 2008).It is in fact the subducting slab-induced cooling that is responsible for this dragging ("viscous blanket" formation; Kincaid and Sacks, 1997, see also Fig. 1), or, in other words, the "lithospherisation" of the asthenosphere at the crust surface that triggers its downwards flow.One may wonder if the interplate geotherm finally governs the of maximum kinematic decoupling.From the surface to z = z dec , temperature within the interplate channel is strongly affected by the friction coefficient, as γ c controls the shearing stress and thus the rate of dissipated energy (Fig. 7a), if strain rate remains roughly controlled by the constant convergence rate.As a consequence, dissipation rate increases from simulation S13 (23.7 µW m −3 , 50 km depth, γ c = 0.034) to simulations S6 and S12 (28.4 µW m −3 , γ c = 0.045), to a maximum in simulation S10 (39 µW m −3 ,γ c = 0.069).This results in a significant temperature difference along the subduction interplate, as illustrated in Fig. 7b, 12 Myr after simulation start (T = 454 • C in simulation S10 at a depth of 50 km and only 346 • C in simulation S13).Hence the brittle behaviour indirectly affects the interplate decoupling depth by modifying the interplate geotherm.The dissipation energy associated with brittle deformation would rather act as an accelerating factor than a decoupling trigger, actually.In simulation S10LG, performed with an extremely low friction coefficient, dissipation rate and temperature along the subduction interface are very low (Fig. 7).The interplate decoupling depth requires a longer subduction duration to reach a stable location, but finally stabilizes at the same depth as the one modelled in simulation S10 (Fig. 4b).I conclude that it is the thermomechanical equilibrium at the down-dip extent of the interplate plane/mantle wedge tip that finally governs the interplate decoupling depth, and not the brittle behaviour occurring at shallow depths.

Conditions favouring z BDT = z dec
The BDT and the limit of kinematic decoupling between the two converging plates occur at the same depth from 10 Myr after subduction initiation in simulations S13 and S13f14, and from 12 Myr in simulations S6, S12 and S10LG (Fig. 4).In all cases the friction coefficient is low (at maximum equal to 0.045 in simulations S6 and S12, Table 2).A (relatively) high friction coefficient entails high shear stresses that exceed ductile stress at shallow depths and minimizes brittle deformation.A very low friction is hence necessary to sustain brittle behaviour at deep levels.In most cases, the transition from brittle to creep deformation is forced by the sudden high temperature increase in the vicinity of the asthenosphere.The BDT depth would thus be controlled by the decoupling depth location, rather than vice versa.Simulation S6 shows that the weak material filling the interplate boundary is mainly brittle all along the converging boundary.Its upper part jumps to the ductile domain exactly at the depth z = z dec (Fig. 8).The basic influence of corner flow observed when z BDT = z dec is also suggested by the comparison of interplate stress between simulations S13 and S13f14 (Fig. 5), where a low asthenospheric strength shallows z dec and yields a BDT shallowing of the same amount.To sum up, the situation where the BDT occurs at the down-dip extent of the interplate plane is conditioned by a low interplate friction coefficient and/or by a low asthenosphere viscosity at the mantle wedge tip.It is however far from being a general case, as initially assumed by Arcay et al. (2007bArcay et al. ( ,a, 2008)).

Interplate dynamics for a young incoming lithosphere
In this section, the initiation of a young and hot lithosphere subduction under a 100-km-thick upper plate is investigated as a function of the rheological parameter set.The model setup is depicted in Fig. 2b.The same rheologies as in the previous section are tested.The influence of the weak layer density is also tested.Before presenting the modelling results, a synthesis from old plate subduction modelling is briefly summarised to discuss which factors are likely to be the most sensitive for the initiation of young lithosphere subduction.

Parameters favouring the interplate kinematic decoupling: key parameters for the initiation of a hot lithosphere subduction
Strain localisation allowing for kinematic decoupling between the two converging plates is carried out from the surface to the BDT depth thanks to the low crustal brittle yield strength.In the ductile part of the interplate plane (for z BDT < z ≤ z dec ), the localisation of deformation is favoured by a high difference in activation energy between the interplate layer and the surrounding mantle.Localisation is also helped by the non-Newtonian strength, decreasing in high strain rate areas, that focuses deformation within highly deforming and low strength layers.Note that the kinematic decoupling at the interplate down-dip extent is promoted if z BDT = z dec , because shear heating contributes to decrease stress.In the vicinity of the asthenosphere, high strain rates along the interplate plane separate in two main branches (Fig. 1b): one at the viscous blanket surface where the asthenosphere is sheared by slab drag, and the second spreading out along the overlapping plate base where the upper part of the corner flow takes place.The kinematic decoupling stops at z = z dec , because strain localisation within the weak crust layer is shifted towards the viscous blanket surface where the hot asthenosphere is very soft.For a young subducting lithosphere, the kinematic decoupling in the interplate ductile part may be hindered by the hot subducting thermal state: On one hand, the downward slab pull that verticalises the convergent motion is reduced, and, on the other hand, the thermal gradient between subducting slab and hot asthenosphere (and thus the viscous drag efficiency) is lessened.To test the first point, different densities of the weak crust layer are tested.The second one might be partly compensated if a low viscosity is added at the mantle wedge tip.

Results: convergence mode of a young subducting lithosphere -influence of rheology
The subduction of a 18Myr-old (∼ 50 km thick) lithosphere under a 100-km-thick lithosphere is performed with the five rheological sets discussed in section 3 (C6, C10, C12, C13, and C14b) with different densities of the weak layer.Addi-  tional tests are performed with a decreased asthenospheric viscosity.The subduction process initiates nicely in simulations S10b, S13b, and S14b.The kinematic decoupling between the two converging plates and the kinematic coupling with the mantle stir up the corner flow and make the subduction process sustainable (Fig. 9a).On the contrary, in simulations S6b, S12b, and S12b2, the subduction process at the interplate down-dip extent slows down and eventually gets jammed along the sublithospheric layer.In these three cases, the interplate plane is locked very soon (in less than 6 Myr after simulation start) while convergence at the downdip extent, initially with a 30 • dipping angle, becomes horizontal (Fig. 9b).As a consequence, the downwards asthenosphere drag at the subducting slab surface stops, and the upper branch of the corner flow necessary to heat and decouple mechanically the slab surface from the overlying plate base is not active.The subducted slab is thermally weakened and finally drips.Fig. 10.BDT versus interplate decoupling depths modelled for young and old incoming lithospheres.The subducting lithosphere ageing at the trench yields simultaneously a BDT deepening and a shallowing of z dec (dashed arrows), except when a hydrous viscosity reduction is modelled in the asthenosphere, resulting in z BDT = z dec at all ages.
In simulation S12b, subduction fails because the weak crust is not able to dig efficiently its own way towards the asthenosphere and accumulates at shallow depths, which jams the interplate, but the crust buoyancy is set to zero in simulation S12b2 without improving the subduction interplate functioning.The interplate jamming is not observed in simulations S6b and S12b2, although performed with a crust density lower than the mantle density (-190 and -80 kg m −3 , respectively, Table 2), but subduction stops as the kinematic coupling between slab and asthenosphere is still never triggered.Interplate geotherms in simulations S6b, S12b, and S12b2 are close to those simulated in simulations S10b, S13b, and S14b.It is hence not the interplate geotherm that governs alone the efficiency of the interplate kinematic decoupling but probably the corresponding strength contrast between mantle and crust.The three simulations showing a jammed subduction are also the three ones with the lowest E a .Supplementary simulations identical to S6b, S12b and S12b2 are performed by including the hydrous strength decrease applied in rheology C13f14 (A 0 in Eq. 2 is replaced by A 0 /14) without improving the interplate decoupling that cannot occur at the interplate down-dip extent.Hence, E a appears as a basic parameter to model a successful initiation of subduction in this setup.A minimum E a may be necessary to compensate for the cold interplate thermal state to sustain the crust weakness with respect to a normal and thence hotter mantle.Above ∼ 67 km depth, high strain rates jump into the basal sublayer of the forearc lithosphere, suggesting that above this depth the forearc mantle becomes the localising layer (Fig. 9b).The slab flattening is likely to further enhance the interplate cooling and strengthening.However, if a high E a might promote the interplate kinematic decoupling in the vicinity of the hot upper sublithospheric layer, note that a too high E a could result in a buoyant crustal layer too soft to remain attached to the slab, if eclogitization is not modelled.

Interplate geometry for a young incoming lithosphere -comparison to an old plate subduction
Figure 10 compares the BDT and the interplate decoupling depths modelled 12 Myr after subduction initiation, for 18-Myr-old and 100-Myr-old subducting lithospheres.Rheologies for which subduction of a hot lithosphere fails are not included.Two additional simulations, S13c and S13f14c, performed with an intermediate lithosphere age of 30 Myr at the trench (Table 2) are displayed.The situation z dec = z BDT is modelled for a young incoming plate only if a hydrous strength reduction is applied at the asthenospheric wedge tip (simulation S13f14c).Generally, the subducting lithosphere ageing deepens a bit the BDT depth for high friction coefficients (by ∼3 km, rheologies C10 and C14b) but significantly if low frictions are imposed (rheology C13, z BDT deepening of 25 km).The deepening effect of the incoming lithosphere ageing is the direct consequence of the interplate geotherm cooling.Note that the weak crust density modifies a bit the interplate dip angle and a lot the deep slab dip, while buoyant crust reduces dip angles (Table 2).The dip angle decrease associated with crust buoyancy cools the interplate domain and, consequently, deepens both z BDT (offset of only 3 km between simulation S13b and S13b2, Fig. 10) and z dec (deepening by 8.5 km).Contrary to z BDT , the decoupling depth shallows with the subducting lithosphere ageing (by ∼ 10 km and 17 km, respectively, for rheologies C13 and C14b, respectively).This may result from the high viscosity of the viscous blanket when the subducting plate is old: low temperatures at the slab surface strengthen the viscous blanket that is then able to mobilise a thicker layer of hot asthenosphere, and emphasises the corner flow effect on the shallow interplate weakening.This effect is small for rheology C10, where z dec is very deep for an old subducting lithosphere and already close to the upper plate thickness, and cannot be much more deepened for a young lithosphere subduction (deepening of ∼3 km between A lith =100 Myr and A lith =18 Myr).Note that the hydrous strength decrease simulated in the asthenospheric wedge shallows z dec also for a young incoming lithosphere (compare simulations S13c and S13f14c).In that case, the significant z dec shallowing rises z dec up to the BDT depth.As a consequence, the lithosphere ageing effect on z dec cancels the BDT deepening for this specific rheology and shallows the BDT, in spite of interplate cooling.Fig. 11.Compressive force through time, applied on the young subducting plate, to sustain a constant convergence rate of 6.5 cm yr −1 , for different rheological models, subducting lithosphere ages, and weak layer densities.

Consequences for subduction initiation and perennity
To test the viability of the subduction system, the compressive force exerted by the kinematic boundary condition to sustain a constant convergence rate of 6.5 cm yr −1 is computed as the deviation of horizontal deviatoric stresses, σ xx , from hydrostatic stresses in an oceanic column of density ρ ref (Christensen, 1992): where z c is the compensation depth (259 km depth).F s > 0 indicates a compressive stress state.Subduction is realistically modelled if F s approximates that of natural tectonic force, that is, not higher than 10 13 N m −1 .This force becomes rapidly excessive when the mechanical decoupling between the two converging plate does not occur, as simulated in simulations S6b and S12b (Fig. 11).Setting to zero the weak layer buoyancy has no effect on the extremely high exerted compression (simulation S12b2).However, a zero density contrast between weak crust and mantle, ρ c , is required to maintain the applied force to subduct to an acceptable level with rheology C13, even if (1) the subduction interplate decoupling properly works and (2) the weak crust friction coefficient is low.The alternative to model a realistic compressive force (but still high, by ∼ 6.7 × 10 12 N/m) is to decrease the asthenosphere strength in the mantle wedge tip (rheology C13f14).Finally, note that the weak layer friction coefficient strongly increases the force resistant to subduction: close to steady-state subduction, F s ∼ 2.8 × 10 12 N/m in simulation S13 (γ c = 0.034, ρ c =0 kg m −3 ) increases to F s ∼ 7.3 × 10 12 N m −1 in simulation S14b (γ c = 0.069, ρ c =0 kg m −3 ).

Friction coefficient along the interplate shear zone
The coefficient of friction along the subduction fault plane is still a matter of debate, and the low to very low values chosen here may be questionable.Numerous studies of convergence argue for friction coefficients higher than 0.2, such as in northern Chile (Delouis et al., 1998), Himalayas (Cattin and Avouac, 2000), Andaman-Sumatra (Cattin et al., 2009), as also indicated by accretionary/non-accretionary wedge dynamics (Lallemand et al., 1994).On the contrary, very low frictional coefficients (0.05 ≤ γ c < 0.09 and even γ c = 0.03) have been invoked to explain low stress state in Cascadia, Kermadec, NE and SW Japan subduction zones (Wang et al., 1995;Wang and Suyehiro, 1999;Wang and He, 1999;von Herzen et al., 2001) or low shear stress estimates from subduction megathrusts (Lamb, 2006).Similarly, Tichelaar and Ruff (1993)  For all subduction zones, the thermal parameter, φ = A lith × v sub × sin θ , is computed using the Submap database (Heuret and Lallemand, 2005;Heuret et al., 2011, http://submap.gm.univ-montp2.fr/),compiling world-wide subduction rates, v sub , subducting plate ages, A lith , and interplate dip angle, θ .(a) The range of maximum depth of seismogenic rupture along the subduction interplate, D z , encountered in world-wide subduction zones compiled by Heuret et al. (2011), is depicted by the blue dotted domain.Green boxes represent depth intervals of tectonic tremors observed in the vicinity of the subduction interplate, possibly located close to z BDT , after (1) Brown et al. (2009), (2) Brown et al. (2010), (3) Ide ( 2012), (4) Peterson and Christensen (2009), (5) Kao et al. (2009), (6) Boyarko and Brudzinski (2010), and (7) Gomberg et al. (2010).(b) Interplate decoupling depth inferred from seismic tomography (purple dots, depth of high lateral contrast in seismic wave propagation velocity and/or in attenuation in the subduction wedge tip), and z dec estimates from heat flow profiles from trench to back-arc (red boxes and diamonds).For Cascadia, Mexico, N Chile, Central Andes, Sunda and Kermadec, the uncertainty in depth mainly comes from the poorly resolved subduction geometry in the vicinity of the interplate down-dip extent.coefficients of effective friction could result from high pore pressure and porosity down to a few tens of kilometer depth, as suggested in NE Japan (Magee and Zoback, 1993) and in the Cascadia subduction zone (Peacock et al., 2011).Moreover, numerical models of dynamic subduction simulating realistic viscoplastic rheology (including sometimes elasticity) reveal that a low strength interplate plane (γ c < 0.1) is required to model realistic convergence of strong lithospheres (Hassani et al., 1997;Hall and Gurnis, 2003;Sobolev and Babeyko, 2005;Tagawa et al., 2007;Gorczyk et al., 2007).Gerya et al. (2008) show that a very high strength contrast between converging plates and the sheared interplate material, favoured by an interplate friction coefficient close to zero, is necessary to model one-sided subduction over a wide range of subducting lithosphere age.Finally, this suggests strong variations of friction properties from a subduction plane to another, depending on the local lithology, stratigraphy and thermal state.One possible way to reconcile high and low friction coefficient estimates might be to test a friction coefficient significantly varying with depth, as a function of temperature, fluid pressure, and metamorphic reactions.

BDT and subduction decoupling depth in nature
A minimum boundary for z BDT can be inferred from the maximum seismogenic depth along the subduction interface.The maximum depth of seismic coupling was initially thought to be restricted to 40 km depth (Ruff and Kanamori, 1983), but more recent studies suggest a broad interval, between 20-55 km (Tichelaar and Ruff, 1993) and up to 35-70 km (Pacheco et al., 1993, based a 19 subduction zones), which was very recently confirmed by a worldwide subduction catalogue analysis (Heuret et al., 2011), showing that the down-dip extent of seismic coupling is well beyond the Solid Earth, 3, 467-488, 2012 www.solid-earth.net/3/467/2012/upper plate Moho depth in 70 % of subduction zones.Recent mega-earthquakes at subduction zones also suggest that seismogenic slip occurs up to 55 km depth (Lay et al., 2012).Figure 12a compares BDT depths modelled in this paper and maximum interplate seismogenic depths in nature, labelled D z .The comparison is based on the subduction thermal parameter, φ, defined by the product of the lithosphere age at trench, A lith , and of the vertical descending slab velocity, v sub × sin θ, where v sub is the down-dip subduction rate normal to the trench and θ is the interplate dip angle.The simulated interplate dip after 12 Myr of subduction is listed in Table 2.The depth of tectonic tremors, whose hypocenters are found to be close to the interplate fault, can also be used as a bound for z BDT , since tectonic tremors are often argued to occur in the transition realm between seismic slip and aseismic creep (e.g.Ide, 2012;Lay et al., 2012, Fig. 12a).
One might expect that hypocenters of tectonic tremors would be located directly beyond D z , but this is not systematically observed, as for instance, in Alaska, where tectonic tremors could be much deeper than D z , or in Mexico and Costa Rica, where tectonic tremors are shallower than D z (Fig. 12a).
Being aware that the relationship between tectonic tremors and the maximum down-dip extent of the seismogenic interplate may not be straightforward and is well beyond the scope of this study, both types of data are here simply used to constrain the interplate modelling.Considering the range of seismogenic constraints as an upper bound for the BDT depth in the present modelling, rheologies C10 and C14b are the best candidates to fit observations for very different subducting plate ages, meaning that the modelled friction coefficient should be close to 0.069.At low subduction thermal parameters (φ < 1000 km), the BDT depth modelled with rheology C13 (and even C13f14) is in agreement with tectonic tremor depths in eastern Alaska, but only rheology C13f14 is in agreement with D z values observed at high φ (> 3000 km).Low friction coefficients and low E a (rheology C12, even possibly C6) also fit D z estimates for cold subductions, but are inadequate to simulate subductions of a young and hot lithosphere.Very low friction coefficients are excluded (γ c ≤ 0.034, rheologies C13 and C10LG) to model the interplate BDT for cold subduction zones.This finally suggests a minimum friction coefficient of 0.045 to simulate z BDT ≤ 65 km if hydrated asthenosphere viscosity decrease is not simulated (rheology C6, Fig. 4b).The occurrence of z BDT = z dec is then unlikely, unless the mantle wedge presents a low viscosity, possibly sustained by the diving slab dehydration and mantle wedge metasomatism, which questions the necessary amount of strength reduction associated with fluid permeation.How can the interplate decoupling depth in this model be constrained?Two main types of geophysical data may be looked at: seismic wave (attenuation) tomography and heat flows profiles.Once again, the comparison between modelling results and observations is based on the subduction thermal parameter.Abers et al. (2006) show that the sharp lateral transition between hot asthenosphere at the mantle tip and the cooled fore-arc nose on top of the cold subducting slab surface is probably imaged by an abrupt lateral increase in seismic wave velocities and/or in seismic wave quality factor, Q, as observed in NE Japan, Cascadia, and Alaska (Zhao et al., 1992(Zhao et al., , 2001;;Stachnik et al., 2004).z dec would then be inferred at a depth of ∼ 45 km at minimum in Cascadia, ∼65 km in Alaska, and possibly around 65 km in NE Japan (Fig. 12b).It is however difficult to interpret definitely the transition from high to low seismic wave velocity in the vicinity of the slab surface in terms of temperature increase close to the asthenospheric mantle wedge, as low seismic velocity anomalies may be related to either partial melting, possibly compatible with high temperatures, or on the contrary to serpentinisation in the fore-arc lithospheric mantle (Bostock and van Decar, 1995;Bostock et al., 2002;Brocher et al., 2003), implying low temperatures (antigorite being stable below ∼ 650 • C, Schmidt and Poli, 1998).Profiles of surface heat flux perpendicular to the trench display in the area located just after the trench the upper plate surface cooling induced by subduction and fore-arc cooling and, moving far away from the trench, the progressive fore-arc re-heating towards the volcanic front with a significant increase on top of the decoupling depth.The location of the maximum heat flux detected just ahead the volcanic arc, combined with the subducting slab geometry imaged by seismic data and/or Wadati-Benioff zone, can then be used to evaluate z dec .Furukawa (1993) shows that the heat flux increase from the trench toward the volcanic arc can be simulated only if z dec is set to 70 km in thermo-kinematic models of the NE Japan subduction zone (Fig. 12b), in agreement with the z dec value inferred from seismic wave tomography.Two additional cold subduction zones (high φ) with a sufficient resolution in surface heat flow can be discussed: Hikurangi (Townend, 1997, and the Global Heat Flow Database from the International Heat Flow commission, http://www.heatflow.und.edu/index2.html,updated 12 January 2011) and Kermadec (von Herzen et al., 2001).In Hikurangi, the maximum heat flow (110 mW m −2 ) is located ∼ 130 km in front of the volcanic front, while in Kermadec it is very close to the volcanic arc (20-35 km, ∼130 mW m −2 ).Using for both the subduction geometry inferred from the EHB catalog (Engdahl et al., 1998, subduction profiles being extracted through the Submap extraction tool, http://submap.gm.univ-montp2.fr/,Heuret et al., 2011), the corresponding z dec should be around 35-40 km in Hikurangi, and encompassed between 50 and 75 km in Kermadec.In the Cascadia subduction zone, the heat flow increase at the volcanic front (∼ 90 mW m −2 , Davies and Lewis, 1984;Lewis et al., 1988Lewis et al., , 1992;;Hyndman and Lewis, 1999), observed ∼200 km away from the trench, would suggest an interplate decoupling depth of ∼ 75 km based on the subduction geometry inferred from Flück et al. (1997).Nevertheless, thermo-kinematic modelling the mantle wedge flow performed by Currie et al. (2004) Hall and Morley, 2004, and the Global Heat Flow Database, see above), around 330 km away from the trench.Assuming the subduction geometry extracted from the EHB catalogue, the interplate decoupling depth would be encompassed between ∼ 50 and 75 km.In Central Andes, the maximum surface heat flow in the fore-arc region (60 mW m −2 ) is located between 200 and 250 km away from the trench (Springer and Forster, 1998), and was modelled by a contact between subducting slab/fore-arc mantle/hot asthenosphere at around 60-80 km depth (Springer, 1999).However, other very high heat flux measured closer to the trench (50 mW m −2 , Hamza and Munoz, 1996), 120-150 km away from the trench, would rather suggest a shallower z dec , encompassed between 40 and 50 km depth.Finally, the maximum heat flow in the forearc area is detected ∼270 km away from the trench in the Mexico subduction zone (110 mW m −2 , Ziagos et al., 1985), and would correspond to a slab-upper plate decoupling between 50 and 60 km depth following the subduction geometry proposed by Currie et al. (2002).The final depth range for z dec estimates is very broad (35-80 km) and does not vary monotonically with φ.Finally, subductions with high thermal parameter (φ ≥ 2000 km) are better reproduced by rheologies shallowing z dec , that is, including moderate to low E a and/or low asthenospheric viscosities (rheologies C12, C6, C13, C13f14, and C14b).Note that none of them is able to model the very shallow z dec observed in New Zealand.This specific case would advocate for a strong viscosity reduction in the hydrated asthenosphere to simulate a z dec shallower than 50 km, i.e. for a strength decreasing factor much higher than 14 (maybe even higher than 50, Arcay et al., 2008).For hot subduction zones, the above rough estimates of decoupling depth also argue for shallower z dec than modelled.This could be obtained either (1) by increasing the activation volume (Fig. 6), but the chosen V a are already rather high; or (2) by lowering E a , but young plate subduction under thick upper lithosphere is then hindered; or (3) by decreasing the asthenospheric wedge viscosity, which eventually appears, at least in this modelling as the most safety tool to model both a sustainable subduction at any subducting plate age and a subduction interplate with realistic characteristics, for both BDT and kinematic decoupling depths.Moreover, the interplate decoupling depth modelled with rheology C13f14 including a hydrous ductile strength reduction is the closest to what is observed in Alaska.Very hot subduction zones (φ < 400 km) are not modelled in this study, but the deepening effect on z dec of young incoming lithosphere would predict a decoupling depth deeper than modelled for φ = 460 km and, as a consequence, much deeper than observed (Fig. 12b).None of the rheological models tested in this paper is able to repro-duce alone the whole range of observed z BDT and z dec .On the whole, heat flow data indicate shallower z dec than modelled in this paper, especially regarding minimum estimates in N Chile, Central Andes, Hikurangi, and Sunda, which might argue once again for a shallowing factor that could be attributed to asthenospheric weakening in presence of fluids.One must remember that the subduction rate is identical in all the presented simulations.Previous studies showed that convergence rates higher than 6.5 cm yr −1 would lead to decoupling depths shallower than simulated here, whereas low convergence rates (1 cm yr −1 ) would deepen z dec to more than 120 km for a 100-Myr-old subducting lithosphere (Arcay et al., 2007b(Arcay et al., , 2008)).I therefore conclude that the z dec modelling proposed here simulates realistic down-dip extents of the subduction interplate plane, and has the basic advantage to model a self-consistent interplate dynamics.

Conclusions
By combining a non-Newtonian viscous rheology and a pseudo-brittle rheology, thermomechanical models are performed to model the long-term equilibrium state of a subduction interplate.The subduction interplate dynamics is shown to be strongly dependent on both ductile and brittle strength parameters.For an old subducting lithosphere, the brittleductile transition depth mainly depends on the friction coefficient and the activation energy of the interplate material; i.e. high friction and low activation energy shallow the BDT.
If the BDT occurs at the kinematic decoupling depth, it is then affected by the depth dependence of the ductile strength and the asthenosphere strength at the mantle wedge tip.The kinematic decoupling depth along the subduction plane is strongly shallowed if the viscosity in the mantle wedge is low, and moderately to sightly shallowed by an increase in activation volume of the ductile strength, and/or a decrease in E a , the difference in activation energy between mantle and interplate material.Deep BDT can be simulated at the depth of interplate kinematics only if very low interplate friction coefficients are modelled, and/or if a decrease in asthenosphere strength, possibly associated with metasomatism, is included.Regarding young lithosphere subduction, a high activation energy contrast between mantle rocks and the interplate medium is necessary to simulate a realistic subduction of a 20-Myr old plate.Finally, both the BDT depth and the decoupling depth depend on the subducting plate age, but are not influenced in the same fashion: cool and old subducting plates deepen the BDT but shallow the interplate decoupling depth.Even if BDT and kinematic decoupling are intrinsically related to different mechanisms of deformation, this work shows that they are able to interact closely.Comparison between modelling results and observations suggests a minimum friction coefficient of 0.045 for the interplate plane, even 0.069 in some cases, to model realistic BDT depths.The modelled z dec is a bit deeper than suggested by geophysical observations.I ultimately conclude that the better way to improve the adjustment to observations may rely on a moderate to strong asthenosphere viscosity reduction in the metasomatised mantle wedge.
Fig. 1.(a) Definition of the interplate decoupling depth, z dec .An outline of high strain rate is schematically represented by the green line.The light-dotted domain depicts the motionless upper lithosphere.The subduction interplate plane, here envisioned as a tangential kinematic discontinuity, is coloured in black.The interplate weak layer, located at the subducting lithosphere surface as an oceanic subducting crust layer, is depicted in grey.The "viscous blanket" refers to the thermal boundary layer formed by asthenospheric cooling at the subducting lithosphere surface(Kincaid and Sacks, 1997).(b) Definition of the brittle-ductile transition depth, z BDT .The brittle-ductile boundary (red line) connects rock elements where the pseudo-brittle strength, ν b , equals the non-Newtonian viscosity, ν v .The medium is modelled as brittle above and ductile below, as sketch in the stress-depth diagram along the interplate plane on the right.The shallowing effect on z BDT of an energy activation decrease, keeping constant the reference viscosity, is sketched in blue.

Fig. 2 .
Fig.2.Boundary conditions for old and young subducting lithospheres and thermal conditions when subduction is initiated.The weak layer geometry imposed at simulation start is depicted in blue.The subducting velocity, v sub , is imposed on the lithosphere in a 832-km-wide and 16-km-deep domain, respectively, counted from the box left-hand side and from the box surface, respectively.Slip is free on the remaining surface at the top, which allows the subduction geometry for evolving freely.One isotherm every 400 • C. The temperature field is constant in the red dashed area, and mimics the lithosphere cooling from formation at the ridge (top left corner) to a chosen lithosphere age, A lith , 400 km away.(a) Boundary conditions and initial state for a converging lithosphere 100 Myr-old at the trench.A lith is set to 90 Myr.(b) Boundary conditions and initial state for a converging lithosphere ∼20-Myr-old at the trench.A lith is set to 10 Myr.

Fig. 3 .
Fig. 3. (a) Simulation S13 (rheology C13, 100-Myr-old subducting lithosphere, Table 2), 16 Myr after subduction initiation.The trench is located at the abscissa x = 0 km.One isotherm (black line) every 200 • C. The interplate sampling line (orange line) joins the points where dissipation energy rates are maximum (see the text for details).(b) Close-up on the subduction interface.Outlines of dissipation energy rate (green lines) are depicted every 10 µW m −3 .(c) Profiles along the interplate sampling line of, from left to right, (1) dissipated energy rate, (2) second invariant of the deviatoric stress tensor, (3) second invariant of the strain rate tensor, and (4) temperature.Note that, in this particular case, the brittle-ductile transition depth, z BDT , is so close to the interplate decoupling depth, z dec , that they are assumed to be equal.
Fig. 4. (a) brittle-ductile transition depth simulated as a function of time elapsed from subduction initiation, for an old subducting lithosphere (simulations S6, S10, S10LG, S12, S13, S13HG, S13f14, and S14 in Table2).The question mark underlines the assumed estimate for z BDT at 12 Myr for rheology C10LG, with an extremely low friction coefficient for the weak layer, which makes the measurement impossible at that time.Minimum and maximum values can nevertheless be measured, as illustrated by the uncertainty bar.The assumed z BDT value at 12 Myr is set to mid-depth in the uncertainty interval, as it also corresponds to the interplate decoupling depth at this time (b).It is indeed likely that z BDT = z dec in this case.(b) Interplate decoupling depth evolution, for an old subducting lithosphere, simulated for the same models as in (a).

Fig. 5 .
Fig. 5. Stress profiles computed along the subduction interface for an old incoming lithosphere, 12 Myr after subduction initiation, for simulations S10 and S14b in (a), and S13 and S13f14 in (b).(c) Interpretative interplate strength profile as a function of the relationship between z BDT and z dec .
Fig.12.brittle-ductile transition depth (a) and interplate decoupling depth (b) modelled in this study (crosses) as a function of the subduction thermal parameter, φ.For all subduction zones, the thermal parameter, φ = A lith × v sub × sin θ , is computed using the Submap database(Heuret and Lallemand, 2005;Heuret et al., 2011, http://submap.gm.univ-montp2.fr/),compiling world-wide subduction rates, v sub , subducting plate ages, A lith , and interplate dip angle, θ .(a) The range of maximum depth of seismogenic rupture along the subduction interplate, D z , encountered in world-wide subduction zones compiled byHeuret et al. (2011), is depicted by the blue dotted domain.Green boxes represent depth intervals of tectonic tremors observed in the vicinity of the subduction interplate, possibly located close to z BDT , after (1)Brown et al. (2009), (2)Brown et al. (2010), (3) Ide (2012), (4)Peterson and Christensen (2009), (5)Kao et al. (2009), (6)Boyarko and Brudzinski (2010), and (7)Gomberg et al. (2010).(b) Interplate decoupling depth inferred from seismic tomography (purple dots, depth of high lateral contrast in seismic wave propagation velocity and/or in attenuation in the subduction wedge tip), and z dec estimates from heat flow profiles from trench to back-arc (red boxes and diamonds).For Cascadia, Mexico, N Chile, Central Andes, Sunda and Kermadec, the uncertainty in depth mainly comes from the poorly resolved subduction geometry in the vicinity of the interplate down-dip extent.

Table 1 .
Parameter names and values.

www.solid-earth.net/3/467/2012/ Solid Earth, 3, 467-488, 2012 D. Arcay: Modelling of subduction interplate dynamics z
fits heat flow profiles when dec is set to 60 km depth.This is deeper than inferred from seismic wave tomography study, showing a low velocity anomaly at the slab surface starting ∼ 45 km depth, and deepening by following the descending slab surface up to 60 km.In the Sunda-Borneo subduction zone, the heat flow increase occurs ∼ 50 km in front of the volcanic front (70 mW m −2 ,