Numerical modelling of lithosphere–asthenosphere interaction in a subduction zone

Earth and Planetary Science Letters 272 (2008) 698–708

Contents lists available at ScienceDirect

Earth and Planetary Science Letters j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / e p s l

Numerical modelling of lithosphere–asthenosphere interaction in a subduction zone M.-A. Bonnardot a, R. Hassani b,⁎, E. Tric a a b

Laboratoire Géosciences Azur, CNRS, Université de Nice Sophia-Antipolis, 250 rue Albert Einstein, 06560 Valbonne, France Laboratoire de Géophysique Interne et Tectonophysique, CNRS, Université de Savoie, 73376 Le Bourget du Lac, France

A R T I C L E

I N F O

Article history: Received 21 May 2007 Received in revised form 3 June 2008 Accepted 5 June 2008 Available online 18 June 2008 Editor: C.P. Jaupart Keywords: 2-D numerical modelling subduction zone solid–fluid coupling stress regime slab dynamics

A B S T R A C T We developed a new 2-D numerical approach to study solid–fluid coupling applied to subduction zones. The lithosphere is characterised by an elastic or elastoplastic behaviour and the asthenosphere by a homogeneous isoviscous fluid. The temperature effects are ignored and viscosity and density are constant in time. The solid and the fluid problem are discretised by the finite elements method (FEM). The same solid code used in Hassani et al. [Hassani, R., Jongmans, D., Chery, J., Study of plate deformation and stress in subduction processes using two-dimensional numerical models, J. Geophys. Res. (1997) 102 17951–17965.] has been used to compute the solution of the solid problem. The Stokes problem is solved by a direct solver with a stabilisation procedure. We used a very simple staggered coupling method where the fluid domain is regularly re-meshing. We observed numerical instabilities when the time step is not sufficiently small, especially when strong coupling between the solid and the fluid occurs. We have tested different configurations where the lithosphere is elastic or elastoplastic and show how the slab geometry, the topography and the stress regime in the plates are affected by the viscous resistance of the mantle. We observed that the asthenosphere viscosity is a fundamental parameter in the subduction process. For subduction with an extensional regime in the upper plate, we observe a linear decrease of the extensional stress as a function of the asthenospheric viscosity. © 2008 Elsevier B.V. All rights reserved.

1. Introduction Subduction systems result from important feedbacks between a strong overriding plate, a strong subducting plate and a viscous convecting mantle. The main difficulty when attempting to model the whole subduction process can be summarised by the following fundamental question, raised by Han and Gurnis (1999): how to take into account the lithosphere and the asthenosphere as two individual mechanical entities, that are thermally linked and interact strongly due to their belonging to the same convective system? Even if the lithosphere is an integrated part of mantle convection, the large viscosity contrast between the lithosphere and the upper mantle associated with the thermal effects induces a mechanical decoupling between these two entities. In accordance with their rheological behaviour, the subduction of the lithospheric plate can be modelled as the interaction between a solid-mechanical subducting lithosphere and a viscous flowing fluid-like mantle. Many studies have been devoted to this problem during these last decades. None of them has proposed a global method (Lux et al., 1979; Hager and O'Connell, 1981; Schmeling and Jacoby, 1981; Davies, 1986; Gurnis and Davies, 1986; Davies, 1988; Tao and O'Connell, 1993; Christensen and Hofmann, 1994; Zhong and Gurnis, 1995; Christensen, 1996; Han and Gurnis, 1999; Schott et al., 2000) with the exception of ⁎ Corresponding author. E-mail address: [email protected] (R. Hassani). 0012-821X/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.epsl.2008.06.009

Morra and Regenauer-Lieb (2006) who proposed a new approach. Aside from the latter, two fundamentally different approaches exist. In one approach, the lithosphere is modelled as a solid medium, where the material behaviour is governed by elastoviscoplastic laws, overlying an inviscid asthenospheric fluid (Poliakov et al., 1993; Giunchi et al., 1994; Hassani et al., 1997; Toth and Gurnis, 1998; Branlund et al., 2001). The asthenosphere acts on the lithosphere through the hydrostatic pressure, which means that the viscous interaction between the subducted plate and the surrounding mantle is neglected. With this approach, we can test different rheological behaviours of the slab and evaluate their effects on the plate deformation and stress during the whole subduction process. In this way, it has been shown that a significant part of the deformation observed in subduction process was governed by elasticity instead of viscosity only, and that a fully viscoelastic rheology was required to comprehensively model slab dynamics (Funiciello et al., 2003a,b). However, a clear limitation of this approach is that the subduction model is developed neglecting any contribution of forces due to mantle flux. In the second approach, the lithosphere is viewed as a viscous fluid governed by the Navier–Stokes equation. A full viscous coupling between the plates and the upper mantle occurs, allowing for the computation of the mantle flow in response to the subducting plate behaviour. Such a method can use viscous and/or viscoplastic rheologies to simulate processes within slabs and does not identify precisely the mechanical behavior of the lithosphere (Stegman et al., 2006; Schellart et al., 2007).

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

Many studies have been published with different techniques to introduce mobile plates into models of mantle convection. In one technique, called kinematic plate models, the plate velocity is imposed as a surface boundary condition (Turcotte and Oxburgh, 1967; Parmentier and Turcotte, 1978; Lux et al., 1979; Davies, 1986; Gurnis and Davies, 1986; Christensen, 1992; Zhong and Gurnis, 1995; Davies, 1995; Insergueix et al., 1997, 1999). Hager and O'Connell (1981) and Davies and Richards (1992) have found that the surface imposed velocity was proportional to η− 1/4 and η2/3, respectively, where η is the mantle viscosity. A critical velocity of the plate exists above which the mantle convection is organised in a single cell as large as the surface boundary. This velocity corresponds to the dynamic balance between forces exerted on the plate and on the fluid mantle. The second technique, called dynamic plate models, uses either lithospheric weak zones or non-Newtonian viscosity with free-slip boundary conditions at the top surface (Schmeling and Jacoby, 1981; Christensen, 1983; Gurnis and Hager, 1988; King and Hager, 1990, 1994; Schmeling and Schott, 1994; Kincaid and Sacks, 1997; Buck and Poliakov, 1998; Moresi and Solomatov, 1998; Trompert and Hansen, 1998; Zhong et al., 1998; Schmeling et al., 1999; Schott et al., 2000; King et al., 2003; Enns et al., 2005; King and Ita, 2005; Stegman et al., 2006; Schellart et al., 2007). Whatever the technique used, the comparison between fully dynamic and kinematic plate formulation shows that temperature structures and slab evolutionary histories are similar when the effective viscosity and surface velocity are nearly identical (Han and Gurnis, 1999). However, even if we know today the importance of plate rheology (Karato et al., 2001; Funiciello et al., 2003a) and of interplate properties (Shemenda, 1994; Hassani et al., 1997) on the slab geometry and strain deformation of the surrounding plate (Heuret and Lallemand, 2005; Lallemand et al., 2005), we do not know very well the behaviour of the slab imposed by the coupling between the lithosphere and the asthenosphere. The 2-D corner flow model proposed by Tovish et al. (1978) shows that the slab dip and the viscous coupling between the slab and the convecting mantle would be strongly correlated. The mantle flow within the corner located between the overriding plate and the down-going slab induces dynamic pressures that tend to uplift the slab. These authors also showed that the slab meets less frictional resistance when descending into a non-Newtonian mantle than into a Newtonian one. Arcay et al. (2005, 2006) modelled the thermal structure of a subduction zone with a viscous mantle including dehydration reactions and a rheology that depends on pressure, temperature, strain rate and water content. They found that a non-Newtonian rheology leads to a warmer slab surface temperature and to a greater thermal erosion at the base of the overlying plate than found in isoviscous corner flow models. The localization of the eroded region depends on the subducting plate thermal state, and its width increases with high convergence rates and low subduction dip angles. Using seismic wave speed and attenuations anomalies several authors (Karato and Spetzler, 1990; Billen and Gurnis, 2001; Billen et al., 2003; Karato, 2003) argued also that the viscosity in the mantle wedge can be up to 1000 times lower than in the surrounding asthenosphere. This decrease should be mainly associated with the slab dehydration, inducing partial melting of the mantle, and could have strong dynamical consequences for the interaction between the lithosphere and the asthenosphere (Hirth and Kohlstedt, 1996; Iwamori, 1998; Schmidt and Poli, 1998; Braun et al., 2000; Billen and Gurnis, 2001; Ulmer, 2001; Billen et al., 2003; Van Keken, 2003; Arcay et al., 2005, 2006; Richard et al., 2006). In the same way, Čadek and Fleitout (2003) have shown that lateral viscosity variations in the top of the upper mantle can have an important effect on the geoid and on dynamic topography. These lateral viscosity variations control the mechanical coupling between the lithospheric plates and the underlying mantle (Ravine and Morgan, 1993; Karpychev and Fleitout, 2000). In the vicinity of the subduction we observe a comparable effect (Royden and Husson, 2006; Husson, 2006).

699

These contributions and others (Van Keken et al., 2002; Kelemen et al., 2003; Conder, 2005; Manea et al., 2005; Peacock et al., 2005) have led, through the study of the wedge rheology or the link between mantle convection and subduction, to new insights regarding the temperature and evolution of subduction zones. But the fluid–solid coupling between the overriding plate, the subducting plate and the upper mantle remains poorly constrained and has not been properly investigated. To our knowledge, only one model exists which takes account the “solid” behaviour of the lithosphere and the “fluid” behaviour of the asthenosphere applied to subduction zone. It is that of Morra and Regenauer-Lieb (2006). They present an interesting and new dynamic approach for solid–fluid coupling applied to a subduction zone. However, in their approach only the slab is modelled and the overriding plate is neglected. The present study aims to improve our quantitative knowledge of the convecting mantle–lithosphere coupling in subduction zone by modelling a 2-D fluid-structure interaction. We are well aware that this 2-D approach over-estimates some parameters like hydrodynamic stresses in the wedge-corner flow or the lifting torque which acts on the slab. But our main goal is to present a new approach using solid–fluid numerical coupling applied to subduction zones and to explore the effects of this interaction on the subduction system. We focus mainly on effects of the asthenosphere and its viscosity on the slab behaviour as well as on the deformation within the overriding plate and on its topography. We present first the numerical approach that we use to describe the evolution of the system considering a homogenous upper-mantle viscosity and a purely elastic or elastoplastic behaviour for the plate. Finally, we present our results and compare them to previous works. 2. Method and modelling set-up 2.1. Assumptions and limitations At the time-scale of interest for (≥10 Ma) and for reasonable viscosity contrasts between the lithosphere and the upper mantle, the slab and the overriding plate can be viewed as a solid body while the upper mantle can be viewed as a viscous fluid. Thus, modelling the evolution of a subduction zone requires the solution of a fluid-structure interaction problem. We make the following assumptions: (1) The temperature effects are ignored. Viscosity and density are then constant in time. (2) For the sake of simplicity we use an elastic or an elastoplastic one-layer oceanic plate and we simplify the mantle rheology using an isoviscous Newtonian fluid. (3) The frictional Coulomb law is used to model the contact between the subducting and the overriding plates. (4) Only the two-dimensional case is considered. This means, amongst others, that toroidal flow of the mantle is ruled out. Accordingly, this study is limited (1) to relatively fast subductions (compared to the thermal time-scale) and (2) to wide slabs for which the lateral flow of the mantle around the slab can be neglected. This second limitation is directly linked to the 2-D numerical approach; it is a strong limitation and does not allow us to study fully the interaction between asthenosphere and lithosphere. Recently, Stegman et al. (2006) have estimated that the toroidal component of the induced flow can be 3–4 times larger than the poloidal component which could provide an explanation for the curvature of trenches. Thus, effects of toroidal flow can change some of the results presented below, especially in case of a slab retreat (Schellart, 2004; Stegman et al., 2006; Capitanio et al., 2007; Schellart et al., 2007). 2.2. Governing equations Let Ωo and Ωs be the physical domains occupied by the overriding plate and the subducting plate, respectively. The governing equations

700

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

of the quasi-static evolution of the lithospheric plates in the solid domain Ωl = Ωo ∪ Ωs are given by 8 < divσ þ ρl g ¼ 0 in Ωl dσ : ¼ Mðσ ; dÞ in Ωl dt

ð1Þ

where σ is the Cauchy stress tensor, g is the vector due to  acceleration : : T gravity, ρl is the lithosphere density, d ¼ 12 ru þ ru is the Eulerian strain rate tensor and u̇ is the velocity vector. d stands for an objective dt time derivative and M for the constitutive law. Although combination of different constitutive laws M (elastic, plastic and viscous rheologies) can be used in Eq. (1) to model the depth-dependent behaviour of the lithospheric plate, we assume in this work that most of the lithospheric strength is retained by its brittle part and as in Hassani et al. (1997), Buiter et al. (2002) and Bonnardot et al. (2008) a one layer plate with an elastic or elastoplastic rheology is used. Thus, the thickness considered in this model corresponds to the effective elastic thickness of the plate and not to its actual one. In addition, the two plates are pushed against each other by use of kinematic boundary conditions on their vertical edges (see Fig. 1), and Signorini contact constraints with Coulomb friction law are used on the contact boundary ∂Ωo ∩ ∂Ωs: 8 : : δun V0; σ n V0; δun σ n ¼ 0; > > : < jσ t jV−μσ n if: δut ¼ 0; : > σ ¼ −μσ δut > if δut ≠0: : : t n jδut j

ð2Þ

where δu̇t and δu̇n are, respectively, the tangential and the normal component (outward unit normal n is considered) of the relative velocity between a point of one plate and its projection onto the other plate, µ is the Coulomb friction coefficient (assumed constant throughout the contact interface). σn and σt are the normal and tangential stresses, respectively. The stationary Stokes problem: (

ηa Δυ −rp þ ρa g ¼ 0 div υ ¼ 0

in

Ωta

in Ωta

ð3Þ

is used to compute, at a given time t, the viscous stress τ = ηa (∇υ + ∇υT), the pressure p and the velocity field υ in the fluid domain Ωta, knowing the viscosity ηa, the density ρa and a set of boundary

conditions (see Fig. 1). The superscript t in the notation Ωta of the fluid domain is used to point out that the solution of a new Stokes problem is needed at each given time station t. The continuity of the vector traction is required (action–reaction principle), on the interface Γla separating the solid (lithosphere) and the fluid (asthenosphere). Because of the viscosity, the fluid and the solid are perfectly stuck to each other, which also implies the continuity of the velocity field on Γla: 

: υ ¼ u on Γ la ðτ−pIÞd n ¼ σd n

on Γ la

ð4Þ

where n is the unit normal on Γla pointing outward the solid domain. 2.3. Numerical method The set of Eqs. (1)–(3)), complemented by the continuity conditions (4)), completely describes the evolution of the system. To tackle such a coupled problem, several strategies exist (e.g. Felippa et al., 2001). Among them, staggered methods (also called partitioned or weakly coupled methods) are the most commonly used because of their flexibility. In these methods the two sub-problems (1) and (3) are solved separately during one discrete time step (rather than solving the total system (1)–(3)–(4)) and a coupling procedure (generally an iterative scheme) is used by enforcing the continuity conditions (4). The method used in this work consists of a staggered method with an explicit procedure which can be summarised by the following flowcharts: (1) Get the initial solid state (u0̇ ,σ0) and the initial fluid state (υ0, τ0, p0) = (0, 0, 0) (2) New time step: tn + 1 = tn + Δt (3) Estimate the traction vector on Γla using the last known fluid stress: f = (τn − pnI) · n (4) Find (u̇n + 1, σn + 1) by solving the solid problem with the traction boundary condition: σn + 1 · n = f on Γla n+1 (5) Define the fluid domain Ωta according to the new position of the interface Γla n+1 (6) Find (υn + 1, τn + 1, pn + 1) by solving the fluid problem in Ωta n+1 n+1 n with the boundary condition: υ = (u − u ) / Δt on Γla (7) Go to (2)

Fig. 1. Schematic representation of the model at the initial time (not to scale) showing the boundary conditions (a) and solid mesh (blue) and the computed fluid mesh (red) after 600 km of convergence.

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

A 2-D solid FEM code (Hassani et al., 1997) is used to compute the solution of the solid problem. This code is based on an explicit dynamic regularisation method (Underwood, 1983; Cundall and Board, 1988) which needs very small time steps δt. The Stokes problem is also discretised by the FEM and solved by a direct solver with a stabilisation procedure allowing the use of linear elements for both the pressure and the velocity field (Quarteroni and Valli, 1997). At each time station tn = nΔt (Δt ≫ δt, n = 1,2,…, N), a global re-meshing of the fluid domain n Ωta is performed (step 5) and the Stokes problem is solved (step 6). This re-meshing is done by the automatic mesh generator BAMG (Hecht, 1997). As a validation of this coupling procedure we present in the Appendix A the test of the sinking cylinder into a viscous liquid. Although this very simple coupling procedure is easy to implement, it is not without drawbacks. Indeed, it suffers from numerical instabilities if the time step is not sufficiently small, especially when strong coupling between the solid and the fluid occurs (high asthenospheric viscosity or, equivalently, high subduction velocity). This procedure can be slightly improved by a fixed point algorithm consisting of a sub-cycling iterations (see Appendix B). However, more complex implicit/semi-implicit schemes or fully coupled methods (e.g. de Hart, 2002; Deparis et al., 2006; Fernández et al., 2007) must

701

be used for strong coupling. In the results presented below, the time step was chosen such that the maximum viscosity that can be used is large enough for our purpose (ηa ≤ 2 · 1020 Pa s). 2.4. Model set-up Fig. 1a displays the initial configuration of the model used. As we don't model a specific subduction zone, geometrical and material parameters of the plate are fixed for all the simulations. The plate is considered to be an old oceanic one with an elastic thickness H = 50 km, a length L = 1700 km, and a density ρl = 3250 kg m− 3. A fault, initially plane, cuts the entire plate at an angle θ = 15°. As already mentioned, a purely elastic or elastoplastic rheology with a Druker– Prager yield criterion is assumed for the plate. The vertical edges are constrained to normal constant velocities: υs = −3 cm yr− 1 for the subducting plate and υo = 0 cm yr− 1 for the overriding one. Elastic parameters are 1011 Pa for Young's modulus and 0.25 for the Poisson ratio. A cohesion of 107 Pa and a friction angle of 30° are used in the case of elastoplastic plates. The upper-mantle domain is limited, vertically, by its bottom boundary at 700 km depth and the base of the lithosphere and,

Fig. 2. Viscous stress intensity (J2(τ)) and velocity field for a viscosity of 5 × 1019 Pa s at different time steps. The simulation was performed with elastic plates. Δu: amount of convergence.

702

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

horizontally, by two vertical edges, sufficiently far (500 km) from the plate boundaries to avoid boundary effects. Although general boundary conditions can be used, we limit this work to the case of a passive mantle i.e. the mantle flow is only induced by the plate motion). The mantle can then flow inside or outside the domain through the vertical edges where only the normal stress is constrained to be equal to the lithostatic pressure, whereas a zero mantle velocity is imposed on the bottom boundary. For the sake of simplicity, only a homogeneous upper-mantle viscosity equal to the asthenospheric one, ηa, is considered (excepted for one of the simulations where a low viscosity zone is added in the subduction wedge). This asthenosphere viscosity is the main parameter of this study but it is necessary to keep in mind that an increase in the viscosity induces, globally, the same result as an increase in the subduction velocity. In what follows we explore the viscosity ηa in the range 1018–2 × 1020 Pa s which is, according to different sources (e.g. Hirth and Kohlstedt, 1996; Kaufmann and Wolf, 1996; Fjeldskaar, 2000; Čadek and Fleitout, 2003), a plausible value bracket for the mantle viscosity below the oceans. The asthenosphere density is fixed to 3200 kg m− 3; the density contrast between lithosphere and asthenosphere is then Δρ = ρl − ρa = 50 kg m− 3. In each numerical experiment a total time span of 20 Ma is considered (600 km of convergence) and about N = 50 coupling time steps Δt are used. Thus, the fluid domain is re-meshed and the Stokes problem is solved whenever the down-going plate has moved 12 km, which is approximately the element size of the solid mesh. During the numerical computation the total number of elements varies approximately between 9500 and 10,700. As an example, the finite elements mesh obtained after 600 km of convergence and for ηa = 5 × 1019 Pa s is shown on Fig. 1b. The time step used for solving the solid problem is δt = 2.5 yr. 3. Results and discussions In the following two subsections the results obtained with elastic plates are presented and discussed. We show how the slab geometry, the topography and the stress regime in the plates are particularly affected by the viscous resistance of the mantle. The third subsection aims to show the effect of an elastoplastic rheology on the plates deformation.

3.1. Slab dip angle A first set of numerical experiments is conducted with a uniform viscosity distribution in the whole upper-mantle. We recall that the rigidity of the plate and the density contrast are fixed for all experiments, although these parameters also have a strong but more obvious effect on the slab dip (see for example (Hassani et al., 1997; Bonnardot et al., 2008)) for the density contrast effect.) Fig. 2 displays the time evolution of the velocity field and the logarithm of the viscous stress intensity (second invariant of τ) in the upper-mantle for a viscosity of ηa = 5 × 1019 Pa s. As expected, the shear stress is maximum at the base of the overriding plate and on the upper surface of the slab, increasing rapidly in the corner. We note also how the flow changes as the slab progresses into the mantle: during the first stages (Fig. 2a and b) the mantle flows round the slab and leaves the spatial domain Ωta by the left edge. As the slab gets closer to the bottom of the upper-mantle domain, a backward flow gradually takes place below the subducting plate (Fig. 2c and d) in an approximately 300 km thick low-stress channel, while the eddy in the arc corner becomes more vigorous preventing the mantle to go outside by the left edge (Fig. 2d). Except for the use of a refined mesh in the vicinity of the corner, we did not make any attempt to take into account the singular aspect of the solution in the corner. However, the computed non lithostatic pressure in the arc corner (Fig. 3b), when the slab reaches the bottom of the upper-mantle, corresponds closely to an 1/r variation (where r is the distance to the trench) predicted by an analytical model assuming a fixed rigid straight slab crossing the whole upper-mantle (Turcotte and Schubert, 1982). The hydrodynamic forces (often called suction forces) induced by the corner flow (Fig. 3a) works against the slab pull force. Fig. 4a compares the geometry of the lower surface of the slab for different values of the viscosity. The deep dip angle of the slab varies from ca 57° for the inviscid case to ca 40° for ηa = 2 × 1020 Pa s (Fig. 4b). Thus, the lesser the asthenosphere viscosity is, the greater the dip angle of the slab is. For a higher viscosity and for a sufficient slab length, the lifting torque produced on the slab by the suction force exceeds the torque produced by the slab pull force resulting in an up-going motion of the

Fig. 3. (a) The velocity field and the non lithostatic pressure in the corner flow after 500 km of convergence and for a mantle viscosity of 5 × 1019 Pa s. (b) The non lithostatic pressure acting on the upper part of the slab when it reaches the upper-mantle bottom. The numerical solution (circles) is close to a theoretical profile in 1/r (solid line) obtained with a rigid straight slab.

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

703

Fig. 4. Results for an elastic slab. (a) Geometry of the underside of the slab after 600 km of convergence and for different values of the viscosity ηa. Corresponding dip angles (b) and topographies (c).

slab until it flattens under the overriding plate. When this occurs, and as the solid and fluid subsystems are considered individually during one time step in the simple staggered algorithm described above, the out-of-balanced torque leads to very large computed slab velocities. This produces, in turn, very high viscous stress in the mantle which forces the slab to move in the opposite direction, and so on. Therefore, this method fails in such case since it produces numerical instabilities, and a fully coupled method is necessary if one wants to study high viscosity effect and flat subductions. The mechanism of the suction force has been proposed (among several others mechanisms) by van Hunen et al. (2004) to explain flat subduction. Using a purely fluid model they showed that the slab becomes flat for a viscosity of 6.5 × 1020 Pa s. We did not try to find the exact lowest viscosity value that produces flat subduction in our model, but numerical instabilities occur for a value barely greater than 2 × 1020 Pa s. These values are close enough despite the differences between these two approaches. However, it is important to keep in mind that the slab width is infinite by assumption. Thus, trenchparallel flow is not allowed and the lifting force due to poloidal flow is probably overestimated in our modeling (see Dvorkin et al., 1993 for a study on narrow width slabs). This suction mechanism might also be invoked to explain the observed difference, in term of slab dip, between ocean–ocean subductions and ocean–continent ones. Indeed, in a statistical analysis, Lallemand et al. (2005) show that slabs dip more steeply (by about 20° on average) beneath oceanic overriding plate than beneath continental ones. Moreover, mean mantle viscosity profiles inferred from geoid and free-air gravity (Čadek and Fleitout, 2003), from interpretation of the post-glacial land emergence (Kaufmann and Wolf, 1996), or from global circulation models (Becker, 2006) show that the sub-oceanic viscosity is in the range 1018–1020 Pa s, one order of magnitude less

than the subcontinental one. We can then expect a net hydrodynamic force in a sub-oceanic mantle wedge of smaller magnitude than in a subcontinental one. We also note on Fig. 4 that the solutions are very close to each other, and to the inviscid case, if the viscosity is less than 1019 Pa s. This means that for a relatively weak mantle, as expected beneath oceans (e.g. Čadek and Fleitout, 2003; Becker, 2006) and in this kind of model for which the geometry is 2-D and the mantle plays only a passive role, the viscous stress acting on lithospheric plates can be neglected, as done by Shemenda (1993) in laboratory modelling or by Hassani et al. (1997) in numerical modelling. However, we should not forget the importance of the third direction which can induce a strong toroidal mantle flow of which one of the geodynamical consequences would be a curvature of the trench (Morra and Regenauer-Lieb, 2006; Schellart et al., 2007). Such an effect would suggest the viscous stress acting on the slab cannot always be neglected. 3.2. Topography, stress regime and back-arc opening As depicted in Fig. 4c, the surface topography of the overriding plate is directly affected by the mantle viscosity: the amplitude of the depression in the fore-arc region decreases when the mantle viscosity is increased. This depression is a direct consequence of the negative buoyancy force of the slab (Δρ N 0). Indeed, no depression is obtained with a zero density contrast, while a negative one induces a positive topography (Shemenda, 1993; Hassani et al., 1997). The formation of this basin is then interpreted as the result of the wheeling motion of the subducting slab during its length increase, which forces the overriding plate to deflect. Owing to the viscous suction this wheeling motion is less large when the viscosity is high resulting in a shallower basin. We note however that the bottom of this basin can reaches

704

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

Fig. 5. Surface topography obtained with a wedge viscosity (ηw) 10 times smaller than the surrounding asthenosphere (dashed line) and comparison with the ones obtained with a uniform asthenosphere. The low viscosity zone (shaded area on the cartoon) extends to 200 km under the bottom of the overriding lithosphere.

unrealistic depths. This is of course a draw-back of our model, which we choose to keep sufficiently simple and in which the slab pull force is time-increasing (as the density contrast is constant) while the mantle viscosity is depth-independent. Although this depression is also observed for the dynamic topography obtained with purely viscous models (Zhong and Gurnis, 1995; Billen and Gurnis, 2001), our results and interpretations disagree with the ones presented by Billen and Gurnis (2001). They showed that the presence of a low viscosity wedge (with a viscosity ten times smaller than the asthenosphere) leads to a basin of a lesser amplitude compared to the case without a low viscosity wedge and that a uniform reduction of the asthenosphere viscosity does not induce the same result. In their model, the negative topography of the surface is created by the hydrodynamic pressure above the slab (suction) which acts to pull the overriding plate down. The basin is then deeper when the viscosity (and then the suction) is high. This is just the opposite of our results. The main difference between their model and ours is, of course, that the elastic strength of elastoplastic plates is taken into account in the present work. The hydrodynamic pressure in the wedge can sustain the slab but is not able to significantly down-warp the overriding plate. Another difference is the presence of the low viscosity zone in the model of Billen and Gurnis (2001) which extends also within the overriding plate. We have tested the effect of a such a zone in the mantle wedge and we have not noticed a marked difference with the case of a uniform reduction of the asthenosphere viscosity: the corresponding solution associated with a wedge viscosity of 1019 Pa s, 10 times smaller than the rest of the asthenosphere, is, in term of slab geometry and fore-arc basin depth, comprised between the solutions given by a uniform viscosity of 1019 Pa s and 1020 Pa s (Fig. 5). Moreover, it is safe to assume that if we also extend the low viscosity zone to the overriding plate by inserting a viscous material within it, the result should be featured by a more pronounced basin since the plate is weakened. Notwithstanding the presence of this fore-arc basin, the slab pull force does not always induce an extensional regime in the overriding plate. It depends on the frictional coupling between the two plates which produces a compressional horizontal stress that can overcome the traction induced by the slab pull force. Fig. 6 shows the averaged   value hsxx i ¼ H1 ∫0H sxx dz of the horizontal deviatoric stress sxx = (2σxx − σyy − σzz) / 3 computed far from the subduction zone and after 600 km of convergence, as a function of the asthenospheric viscosity and friction coefficient. For a zero Coulomb friction coefficient, bsxxN is positive whatever is the viscosity (in the range [0,2 × 1020 Pa s]), denoting an extensional stress in the overriding plate. Note, however, the (linear) decrease of the extensional stress with ηa (Fig. 6). Consequently, if a zone of weakness (or a localisation mechanism) exists in the overriding plate enabling a back-arc basin to open (Hassani et al., 1997; Buiter et al., 2002), the opening would be less in case of a strong viscous asthenosphere, and in fact inhibited if there is a frictional resistance (even low) between the two plates. According to

these results, we conclude that a positive density contrast Δρ, a low friction coefficient µ and a low asthenospheric viscosity ηa are factors that promote extensional back-arc regime. 3.3. Elastoplastic rheology As in Hassani et al. (1997) we mimic the non linear and pressuredependent behaviour of the lithosphere by using an elastoplastic rheology with a Drucker–Prager criterion assuming a cohesion c = 10 MPa and a friction angle ϕ = 30° for the two plates: f ðσ Þ ¼ J2 ðσ Þ þ α ð/ÞðI1 ðσ Þ−p0 ÞV0

qffiffi where I1 ðσ Þ ¼ 13 traceðσ Þ, J2 ðσ Þ ¼ 32jjσ −I1 ðσ ÞIjj, α = 6sinϕ / (3 − sinϕ) and p0 = c / tan(ϕ). All others rheological and geometrical parameters are as in Section 3.1. Fig. 7 shows the dramatic difference between the result obtained with elastoplastic plates and that obtained with elastic ones. For an inviscid asthenosphere (ηa = 0) and for elastoplastic plates, the maximum amount of convergence we can achieve (without re-meshing the plates domain) is slightly greater than 400 km. Indeed, due to plastic dissipation the down-going plate does not unbend like in the purely elastic case and this results in a higher slab dip and then a more deformed overriding plate. Moreover, the strain pattern (Fig. 7) shows a highly localised strain in the fore-arc zone and in the top slab region, which can be viewed as the zones where the two plates may break off. The unrealistic topography (see also Fig. 8c) we obtain with a low asthenosphere viscosity indicates that one of the two plates should probably break off before the overriding plate reaches this exaggerated topography.

  Fig. 6. Averaged horizontal deviatoric stress hsxx i ¼ H1 ∫0H sxx dz within the overriding plate, computed far from the subduction zone and after 600 km of convergence, versus asthenosphere viscosity ηa, for three values of the friction coefficient µ. Positive values of bsxxN indicate extension. The set of symbols corresponds to the different experiments and dotted lines are linear fits. Extensional stress decreases with both ηa and µ. For µ = 0.1, for example, the stress regime becomes compressional from a viscosity of 5 × 1019 Pa s, while for this viscosity an extensional stress greater than 80 MPa is obtained for µ = 0.05.

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

705

Fig. 7. Second invariant J2(e)of the total deviatoric strain e(t) = ∫0t d(τ)dτ obtained after 400 km of convergence (t = 13.3 Ma) with an inviscid asthenosphere and for elastoplastic plates (a) and elastic ones (b).

We again conduct a set of five experiments by changing the asthenosphere viscosity ηa in the range [0,2 × 1020] Pa s. The results are displayed in Fig. 8. The deep slab dip angle is close to 70° for ηa = 0 Pa s

while for ηa = 2 × 1020 Pa s is 40°, approximately. However, the dip angle variations with depth (Fig. 8b) are not as smooth as they are with an elastic slab (Fig. 4b).

Fig. 8. Results for an elastoplastic slab. As in Fig. 5, slab geometry (a), dip angle (b) and surface topography (c) are displayed for five values of the asthenospheric viscosity ηa.

706

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

4. Conclusions We have proposed a new numerical approach for solid–fluid mechanical coupling applied to subduction zones. In this approach, the lithosphere is characterised by an elastic or elastoplastic behaviour, and the mantle by a homogeneous isoviscous fluid. The thermal effects are not taken into account. The method used in this work consists of a staggered method with an explicit procedure between the “solid” solver and “fluid” solver. A re-meshing of the domains is done by the automatic mesh generator BAMG. This approach requires small time steps especially when a strong coupling between the solid and the fluid occurs, which may favour numerical instabilities. The procedure has been slightly improved by a fixed algorithm consisting of a sub-cycling iteration. However, a more robust coupling method (but not so easy to implement) should be used if one wants to study faster subduction or slab interaction with a more viscous material (interaction with the lower mantle for example). Two other promising approaches could be more suitable for this problem, especially in a 3-D context. The first one is the so called Fictitious Domain Method (FDM) where the fluid problem is solved on a fixed fictitious domain which includes both the real fluid domain and the solid one. Lagrange multipliers are used to enforce the continuity of the velocity and traction fields on the solid–fluid interface. Such a method avoids the need for re-meshing (which is a difficult and time-consuming task in 3-D) but coupling strategies are still required to overcome numerical instabilities. The second one is the method described in Morra and Regenauer-Lieb (2006) where a Boundary Elements Method (for the fluid problem) is coupled to a Finite Elements Method (for the solid problem) through the calculation of a drag tensor. The authors note that this method allows for a more stable coupling than staggered methods. Currently, we did not try this interesting approach and only the FDM technique was tested in 2-D. Results obtained with the FDM on some benchmarks are comparable with those obtained by the method presented in this work; we think however that the main advantage of the FDM will be better highlighted in 3-D. We have tested feedbacks between strong overriding and subducting plates and the asthenosphere as a function of the mantle viscosity and/or the rheological behaviour of the lithosphere. Whatever the rheological behaviour of the lithosphere, we observed that the lower the asthenosphere viscosity is, the greater the dip angle of the slab is. These results are associated with a more or less important variations of the topographic surface with a strong localisation of the strain in the fore-arc zone and in the top slab region for an elastoplastic behaviour. A significant result of our study concerns the stress regime in the overriding plate. We observe that the extensional stress in the overriding plate is a linear decreasing function of the asthenospheric viscosity for a given friction coefficient. These results confirm that the overriding plate deformation, the slab deformation and the slab geometry are mainly linked to the asthenospheric viscosity that constitutes the main parameter of the coupling between the slab and the overriding plate. We are well aware that an isoviscous fluid is not representative of a realistic mantle, but it is a first step in the development of a more complete numerical code in which the thermal effects and a non-Newtonian fluid will be included. In its current version, our model also makes it possible to carry out real comparisons with laboratory models like those developed by A. Chemenda, C. Faccenna, F. Funiciello, J. Martinod, V. Regard, and others (Shemenda, 1993; Faccenna et al., 1996; Funiciello et al., 2003a; Schellart, 2005; Regard et al., 2006; Heuret et al., 2007), which will constitute a next contribution. Acknowledgments This research was supported by the CNRS-INSU DyETI program Dynamics of Subduction. We are grateful to Guust Nolet for a detailed

reading of the manuscript. Thanks to Gabrielle Morra and Wouter Schellart, the two reviewers, for their useful and constructive comments. Appendix A . The sinking cylinder As a simple validation test we use the example of the falling cylinder in a viscous fluid channel. Let R be the radius of the cylinder, ρc its density, 2L the channel width, ρf the density of the fluid and η its viscosity. According to (Wang and Liu, 2004), if the cylinder is rigid and the fluid channel infinitely long, the terminal velocity reached by the cylinder is given by: Vlim ¼

 ðρc −ρf ÞgR2  1:72442 −1:73024 −ln−0:9157 4η

with ξ = R / L, while, of course, the total drag force (per unit span) tends to the effective weight of the cylinder: Flim = (ρc − ρf)gπR2. The test is performed with thepfollowing parameters: ρc = 3000 kg m− 3, ffiffiffi ρf = 1000 kg m− 3, R ¼ 2mm, L = 4R and η = 1 Pa s. We can solve the dynamical problem by adding to the right hand side of the equation of motion in (1) the inertial volume force ρcü. As in Wang and Liu (2004) we examine two cases: the case of a rigid cylinder and the case of a very soft one. To simulate a rigid behaviour all nodes of the solid mesh are forced to have the same trajectories, i.e. we simply solve the equation of motion for only one degree of freedom. The very deformable cylinder is simulated by taking a very low Young's modulus: E = 1 kPa. In the latter no exact solution exists to the best of our knowledge. The only thing that we know is that the terminal velocity should be greater than the terminal velocity of the rigid cylinder. The Fig. 9a shows that the computed velocity compares well with the theoretical terminal velocity in the rigid case while, as expected for

Fig. 9. (a) Time evolution of the cylinder velocity for the rigid case (continuous line) and for the deformable case (dotted line). The dashed line corresponds to the theoretical terminal velocity for the case of a rigid cylinder in infinitely long channel. (b) Distribution of the viscous stress intensity (in Pa) near the cylinder and the corresponding velocity field at a given time and for a deformable cylinder.

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

the deformable case, the velocity is increased by the flexibility of the cylinder. The velocity field and the distribution of the deviatoric viscous stress intensity at a given time and in the vicinity of the deformable cylinder are depicted on the Fig. 9b. Appendix B. Numerical instabilities example We have mentioned that, unfortunately, in some cases the explicit coupling procedure suffers from numerical instabilities. The procedure can be slightly improved by using a fixed point algorithm consisting of sub-cycling iterations, repeating instructions (3)–(6) until no changes hold. The new algorithm can be summarised as following: 0 ̇

0

0

(1') Get the initial solid state (u , σ ) and the initial fluid state (υ , τ0, p0) = (0, 0, 0) (2') New time step: tn + 1 = tn + Δt. Set k = 0 and (τn + 1,0, pn + 1,0) = (τn, pn) (3') Estimate the traction vector on Γla using the last known fluid stress: f = (τn + 1,k − pn + 1,kI) · n (4') Find (u̇n + 1,k, σn + 1,k) by solving the solid problem with the traction boundary condition: σn + 1,k · n = f on Γla n + 1,k (5') Define the fluid domain Ωta according to the new position of the interface Γla n + 1,k (6') Find (υn + 1,k, τn + 1,k, pn + 1,k) by solving the fluid problem in Ωta n+1 n + 1,k n with the boundary condition: υ = (u − u ) / Δt on Γla. Set k := k + 1 and go to (3′) if necessary In a less time-consuming variant, the fluid domain is defined only n + 1,0 during the sub-cycling once a time step by taking it to be Ωta iterations on the time step. This avoids the re-meshing of this domain at each iteration and the computation/factorization of the corresponding finite element matrix. We present an example where numerical instabilities take place and show the improvement made by the fixed point algorithm. The numerical experiment consists of an elastic plate of length l = 3 m and

Fig. 10. A simulation used to exhibit numerical instabilities. An elastic plate embedded at its top end and immersed in a viscous fluid is bent by the flow where the inflow velocity is gradually increased until a value of 60 m s− 1 is reached. (a) Fluid pressure and velocity field at the end of the simulation. (b) Time evolution of the total viscous force that acts on the plate computed with the fully explicit procedure (solid line) and with the modified procedure (dashed line).

707

thickness h = 0.2 m immersed in a fluid channel of width L = 6 m. The fluid, of viscosity η = 105 Pa s, flows into the channel with a prescribed horizontal velocity which increases linearly from 0 to 60 m/s between 0 and 0.5 s and then remaining constant. Consequently, the plate is increasingly bent until the fluid velocity reaches its maximal value. The fluid pressure and the velocity field at the end of the simulation are displayed on the Fig. 10a. The two components of the total drag force (per unit length) exerted by the fluid flow on the surface plate Γ, i.e.: F ¼ ∫Γ σndS are shown on Fig. 10b for the two procedures. With the initial procedure and as the plate is supposed to reach a stationary state beyond t = 0.5 s numerical instabilities develop, inducing an oscillating motion of the plate. In this example the modified procedure is able to overcome this problem; however, it does not seems to be a robust method since increasing the viscosity or the inflow velocity will produce again numerical instabilities. References Arcay, D., Tric, E., Doin, M.-P., 2005. Numerical simulations of subduction zones. Effect of slab dehydration on the mantle wedge dynamics. Phys. Earth Planet. Inter. 149, 133–153. Arcay, D., Tric, E., Doin, M.-P., Bousquet, P., de Capitani, C., 2006. Overriding plate thinning in subduction zones: localized convection induced by slab dehydration. Geochem. Geophys. Geosyst. 7, 1–26. Becker, T.W., 2006. On the effect of temperature and strain-rate dependent viscosity on global mantle flow, net rotation, and plate-driving forces. Geophys. J. Int. 167, 943–957. Billen, M.I., Gurnis, M., 2001. A low viscosity wedge in subduction zones. Earth Planet. Sci. Lett. 193, 227–236. Billen, M.I., Gurnis, M., Simons, M., 2003. Multiscale dynamics of the Tonga–Kermadec subduction zones. Geophys. J. Int. 153, 359–388. Bonnardot, M.A., Hassani, R., Tric, E., Ruellan, E., Régnier, M., 2008. Effect of margin curvature on plate deformation in 3-D numerical model of subduction zones. Geophys. J. Int. 173, 1084–1094. Branlund, J.M., Regenauer-Lieb, K., Yuen, D.A., 2001. Weak zone formation for initiating subduction from thermo-mechanical feedback of low-temperature plasticity. Earth Planet. Sci. Lett. 190, 237–250. Braun, M., Hirth, G., Parmentier, E., 2000. The effect of deep damp melting on mantle flow and melt generation beneath mid-ocean ridge. Earth Planet. Sci. Lett. 176, 339–356. Buck, W.R., Poliakov, A.N.B., 1998. Abyssal hills formed by stretching oceanic lithosphere. Nature 392, 272–275. Buiter, S.J.H., Govers, R., Wortel, M.J.R., 2002. Two dimensional simulations of surface deformation caused by slab detachment. Tectonophysics 354, 195–210. Čadek, O., Fleitout, L., 2003. Effect of lateral viscosity variations in the top 300 km on the geoid and dynamic topography. Geophys. J. Int. 152, 566–580. Capitanio, F.A., Morra, G., Goes, S., 2007. Dynamic models of downgoing plate-buoyancy driven subduction: subduction motions and energy dissipation. Earth Planet. Sci. Lett. 262, 284297. Christensen, U.R., 1983. Convection in a variable-viscosity fluid: Newtonian versus power-law rheology. Earth Planet. Sci. Lett. 64, 153–162. Christensen, U.R., 1992. An Eulerian technique for thermomechanical modeling. J. Geophys. Res. 97, 2015–2036. Christensen, U.R., 1996. The influence of trench migration on slab penetration into the lower mantle. Earth Planet. Sci. Lett. 140, 27–39. Christensen, U.R., Hofmann, A.W., 1994. Segregation of subducted oceanic crust in the convecting mantle. J. Geophys. Res. 99, 19867–19884. Conder, J.A., 2005. A case for hot slab temperatures in numerical viscous flow models of subduction zones with an improved fault zone parametrization. Phys. Earth Planet. Inter. 149, 155–164. Cundall, P.A., Board, M., 1988. A microcomputer program for modeling large strain plasticity problems. In: swoboda, G. (Ed.), Numerical Methods in Geomechanics. Balkena, Rotterdam. Davies, G.F., 1986. Mantle convection under simulated plates: effects of heating modes and ridge and trench migration, and implications for the core–mantle boundary, bathymetry, the geoid and Benioff zones. Geophys. J. R. Astron. Soc. 84, 153–183. Davies, G.F., 1988. Role of the lithosphere in mantle convection. J. Geophys. Res. 93, 10451–10466. Davies, G.F., 1995. Penetration of plates and plumes through the mantle transition zone. Earth Planet. Sci. Lett. 133, 507–516. Davies, G.F., Richards, M.A., 1992. Mantle convection. J. Geol. 100, 151–206. de Hart, J., Fluid-structure interaction in the aortic heart valve. A three-dimensional computational analysis, PhD thesis, University of Eindhoven, The Netherlands, 2002. Deparis, S., Discacciati, M., Fourestey, G., Quarteroni, A., 2006. Fluid-structure algorithms based on Steklov–Poincaré operators. Comput. Methods Appl. Mech. Eng. 195, 5797–5812.

708

M.-A. Bonnardot et al. / Earth and Planetary Science Letters 272 (2008) 698–708

Dvorkin, J., Nur, A., Mvako, G., Ben-Avraham, Z., 1993. Narrow subducting slabs and the origin of backarc basins. Tectonophysics 227, 63–79. Enns, A., Becker, T.W., Schmelling, H., 2005. The dynamics of subduction and trench migration for viscosity stratification. Geophys. J. Int. 160, 761–775. Faccenna, C., Davy, P., Brun, J.-P., Funiciello, R., Giardini, D., Mattei, M., Naplas, T., 1996. The dynamic of backarc basins: an experimental approach to the opening of the Thyrrhenian Sea. Geophys. J. Int. 126, 781–795. Felippa, C.A., Park, K.C., Farhat, C., 2001. Partitioned analysis of coupled mechanical systems. Comput. Methods Appl. Mech. Eng. 190, 3247–3270. Fernández, M.A., Gerbeau, J.-F., Grandmont, C., 2007. A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid. Int. J. Numer. Methods Eng. 69 (4), 794–821. Fjeldskaar, W., 2000. What about the asthenosphere viscosity? Comment on ‘Sea-level change, glacial rebound and mantle viscosity for northern Europe’ by K. Lambeck, C. Smither and P. Johnston. Geophys. J. Int. 142, 277–278. Funiciello, F., Morra, G., Regenauer-Lieb, K., Giardini, D., 2003a. Dynamics of retreating slabs: 1. Insights from two-dimensional numerical experiments. J. Geophys. Res. 108. Funiciello, F., Morra, G., Regenauer-Lieb, K., Giardini, D., 2003b. Dynamics of retreating slabs: 2. Insights from three-dimensional laboratory experiments. J. Geophys. Res. 108. Giunchi, C., Gasperini, P., Sabadini, R., D'Agostino, G., 1994. The role of subduction on the horizontal motions in the Thyrrhenian basin: a numerical model. Geophys. Res. Lett. 21, 529–532. Gurnis, M., Davies, G.F., 1986. Mixing in numerical models of mantle convection incorporating plate kinematics. J. Geophys. Res. 91, 6375–6395. Gurnis, M., Hager, B.H., 1988. Controls on the structure of subducted slabs and the viscosity of the lower mantle. Nature 335, 317–321. Hager, B.H., O'Connell, R.J., 1981. A simple global model of plate dynamics and mantle convection. J. Geophys. Res. 86, 4843–4867. Han, L., Gurnis, M., 1999. How valid are dynamical models of subduction and convection when plate motions are prescribed? Phys. Earth Planet. Inter. 110, 235–246. Hassani, R., Jongmans, D., Chery, J., 1997. Study of plate deformation and stress in subduction processes using two-dimensional numerical models. J. Geophys. Res 102, 17951–17965. Hecht, F., 1997 BAMG: bidimensional anisotropic mesh generator, http://pauillac.inria. fr/cdrom/prog/unix/bamg/fra.htm, INRIA. Heuret, A., Lallemand, S., 2005. Plate motions, slab dynamics and back-arc deformation. Phys. Earth Planet. Inter. 149, 31–51. Heuret, A., Funiciello, F., Faccenna, C., Lallemand, S., 2007. Plate kinematics, slab shape and back-arc stress: a comparison between laboratory models and current subduction zones. Earth Planet. Sci. Lett. 256, 473–483. Hirth, G., Kohlstedt, D.L., 1996. Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere. Earth Planet. Sci. Lett. 144, 93–108. Husson, L., 2006. Dynamic topography above retreating subduction zones. Geology 34, 741–744. Insergueix, D., Dupeyrat, L., Menvielle, M., Tric, E., 1997. Dynamical and thermal structure of subduction zone: influence of slab geometry on the convective state of the Earth's upper mantle. Preliminary results. Phys. Earth Planet. Inter. 99, 231–249. Insergueix, D., Dupeyrat, L., Tric, E., Menvielle, M., 1999. Influence of plate kinematics, subduction geometry and convection intensity on Earth's upper mantle dynamics in the vicinity of a subduction zone. Geophys. J. Int. 138, 275–284. Iwamori, H., 1998. Transportation of H2O and melting in subducting zones. Earth Planet. Sci. Lett. 160, 65–80. Karato, S.-I., 2003. Mapping water content in the upper mantle. In: Eiler, J. (Ed.), Subduction Factory. Geophys. Monogr. Ser. AGU, Washington, DC. Karato, S.-I., Spetzler, H., 1990. Defect microdynamics and physical mechanisms of seismic wave attenuation and velocity dispersion in the Earth's mantle. Rev. Geophys. 28, 399–421. Karato, S.-I., Riedel, M.R., Yuen, D.A., 2001. Rheological structure and deformation of subducted slabs in the mantle transition zone: implications for mantle circulation and deep earthquakes. Phys. Earth Planet. Inter. 127, 83–108. Karpychev, M., Fleitout, L., 2000. Simple consideration on forces driving plate motions and on the plate-tectonic contribution to the long-wavelength geoid. Geophys. J. Int. 127, 268–282. Kaufmann, G., Wolf, D., 1996. Deglacial land emergence and lateral upper-mantle heterogeneity in the Svalbard Archipelago-II. Extended results for high-resolution load models. Geophys. J. Int. 127, 125–140. Kelemen, P.B., Rilling, J.L., Parmentier, E.M., Mehl, L., Hacker, B.R., 2003. Thermal structure due to solid-state flow in the mantle wedge beneath arcs. In: Eiler, J. (Ed.), Inside the Subduction Factory. Geophys. Monogr., Ser., vol. 138. AGU, Washington, D.C., pp. 293–311. Kincaid, C., Sacks, I.S., 1997. Thermal and dynamical evolution of the upper mantle in subduction zones. J. Geophys. Res. 102, 12295–12315. King, S.D., Hager, B.H., 1990. The relationship between plate velocity and trench viscosity in Newtonian and power-law subduction calculations. Geophys. Res. Lett. 17, 2409–2412. King, S.D., Hager, B.H., 1994. Subducted slabs and the geoid. 1. Numerical experiments with temperature-dependent viscosity. J. Geophys. Res. 99, 19843–19852. King, S.D., Ita, J.J., 2005. The effect of slab rheology on mass transport across a phase transition boundary. J. Geophys. Res. 100, 20211–20222. King, S.D., Gable, C.W., Weinstein, S.A., 2003. Models of convection-driven tectonic plates: a comparison of method and results. Geophys. J. Int. 109, 481–487. Lallemand, S., Heuret, A., Boutelier, D., 2005. On the relationships between slab dip, back-arc stress, upper plate absolute motion, and crustal nature in subduction zones. Geochem., Geophys., Geosyst. 6.

Lux, R.A., Davies, G.F., Thomas, J.H., 1979. Moving lithospheric plates and mantle convection. Geophys. J. R. Astron. Soc. 57, 209–228. Manea, V.C., Manea, M., Kostoglodov, V., Sewell, G., 2005. Thermo-mechanical model of the mantle wedge in Central Mexican subduction zone and a blob tracing approach for the magma transport. Phys. Earth Planet. Inter. 149, 165–186. Moresi, L.N., Solomatov, V.S., 1998. Numerical investigation of 2D convection with extremely large viscosity variations. Phys. Fluids 7, 2154–2162. Morra, G., Regenauer-Lieb, K., 2006. A coupled solid–fluid method for modeling subduction. Philos. Mag. 86, 3307–3323. Parmentier, E.M., Turcotte, D.L., 1978. Two-dimensional mantle flow beneath a rigid accreting lithosphere. Phys. Earth Planet. Inter. 17, 281–289. Peacock, S.M., Van Keken, P.E., Holloway, S.D., Hacker, B.R., Abers, G.A., Fergason, R.L., 2005. Thermal structure of the Costa Rica–Nicaragua subduction zone. Phys. Earth Planet. Inter. 149, 187–200. Poliakov, A.N., Cundall, P.A., Y Podladchikov, Y., Lyakhovsky, V.A., 1993. An explicit inertial method for the simulation of the viscoelastic flow: an evaluation of elastic effects on diapiric flow in two- and three-layers models. In: Stone, D.B., Runcorn, S.K. (Eds.), Flow and Creep in the Solar System: Observations, Modelling and Theory. Kluwer Academic Publishers, pp. 175–195. Quarteroni, A., Valli, A., 1997. Numerical Approximation of Partial Differential Equations, 2nd ed. Springer. Ravine, M.A., Morgan, P., 1993. Geoid effects of lateral viscosity variations near the top of the mantle: a 2-D model. Earth Planet. Sci. Lett. 119, 617–625. Regard, V., Faccenna, C., Martinod, J., Bellier, O., 2006. Slab pull and indentation tectonics: insights from 3D laboratory experiments. Phys. Earth Planet. Inter. 149, 99–113. Richard, G., Bercovici, D., Karato, S.-I., 2006. Slab dehydration in the Earth's mantle transition zone. Earth Planet. Sci. Lett. 251, 156–167. Royden, L.H., Husson, L., 2006. Trench motion, slab geometry and viscous stresses in subduction systems. J. Geophys. Int. 167, 881–905. Schellart, W.P., 2004. Kinematics of subduction and subduction-induced flow in the upper mantle. J. Geophys. Res. 109, B07401. doi:10.1029/2004JB002970. Schellart, W.P., 2005. Influence of the subducting plate velocity on the geometry of the slab and migration of the subduction hinge. Earth Planet. Sci. Lett. 231, 197–219. Schellart, W.P., Freeman, J., Stegman, D.R., Moresi, L., May, D., 2007. Evolution and diversity of subduction zones controlled by slab width. Nature 446, 308–311. Schmeling, H., Jacoby, W.R., 1981. On modelling the lithosphere in mantle convection with non-linear rheology. J. Geophys. 50, 89–100. Schmeling, H., Schott, B., 1994. On the dynamics of lithospheric roots DFT workshop on modelling of orogenic processes at lithospheric scale, March 21–23, Berlin. Terra Nostra 93, 139–141. Schmeling, H., Monz, R., Rubie, D.C., 1999. The influence of olivine metastability on the dynamics of subduction. Earth Planet. Sci. Lett. 165, 55–66. Schmidt, M., Poli, S., 1998. Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation. Earth Planet. Sci. Lett. 163, 361–379. Schott, B., Yuen, D.A., Schmeling, H., 2000. The diversity of tectonics from fluiddynamical modeling of the lithosphere–mantle system. Tectonophysics 322, 35–51. Shemenda, A.I., 1993. Subduction of the lithosphere and back-arc dynamics: insights from physical modeling. J. Geophys. Res. 98, 16167–16185. Shemenda, A.I., 1994. Subduction: insights from physical modelling. Ser. Modern Approaches in Geophysics, Netherlands. Kluwer Academic Publishers, p. 215. Stegman, D.R., Freeman, J., Schellart, W.P., Moresi, L., May, D., 2006. Influence of trench width on subduction hinge retreat rates in 3-D models of slab rollback. Geochem., Geophys., Geosyst. 7, Q03012. doi:10.1029/2005GC001056. Tao, W.C., O'Connell, R.J., 1993. Deformation of a weak slab and variation of seismicity with depth. Nature 361, 626–628. Toth, J., Gurnis, M., 1998. Dynamics of subduction initiation at preexisting fault zones. J. Geophys. Res. 103, 18053–18067. Tovish, A., Schubert, G., Luyendyk, B.P., 1978. Mantle flow pressure and the angle of subduction: non-Newtonian corner flows. Geophys. J. Int. 83, 5892–5898. Trompert, R., Hansen, U., 1998. Mantle convection simulations with rheologies that generate plate-like behavior. Nature 395, 686–689. Turcotte, D.L., Oxburgh, E.R., 1967. Finite amplitude convection cells and continental drift. J. Fluid Mech. 28, 29–42. Turcotte, D.L., Schubert, G., 1982. Geodynamics. Application of Continuum Physics to Geological Problems. John Wiley, New York. 450 pp. Ulmer, P., 2001. Partial melting in the mantle wedge — the role of H2O in the genesis of mantle-derived ‘arc-related’ magmas. Phys. Earth Planet. Inter. 127, 215–232. Underwood, P., 1983. Dynamic relaxation. In: Belytschko, T., Hughes, T.J.R. (Eds.), Computational Methods for Transient Analysis. Elsevier Sci., New York. van Hunen, J., van den Berg, A., Vlaar, N.J., 2004. Various mechanisms to induce presentday shallow flat subduction and implications for the younger Earth: a numerical parameter study. Phys. Earth Planet. Inter. 146, 179–194. Van Keken, P.E., 2003. The structure and dynamics of the mantle wedge. Earth Planet. Sci. Lett. 215, 323–338. Van Keken, P., Kiefer, B., Peacock, S., 2002. High-resolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle. Geochem., Geophys., Geosyst. 3. Wang, X., Liu, W.K., 2004. Extended immersed boundary method using FEM and RKPM. Comput. Methods Appl. Mech. Eng. 193, 1305–1321. Zhong, S., Gurnis, M., 1995. Mantle convection with plates and mobile, faulted plate margins. Science 267, 838–843. Zhong, S., Gurnis, M., Moresi, L., 1998. Role of faults, nonlinear rheology, and viscosity structure in generating plates from instantaneous mantle flow models. J. Geophys. Res. 103, 15255–15268.