<figure><img src="image_1.png" /></figure>

HAL
open science

# Dependency of slab geometry on absolute velocities and conditions for cyclicity: insights from numerical modelling,

M. Gerbault, E. Tric, R. Hassani, G. Gibert

► To cite this version:

M. Gerbault, E. Tric, R. Hassani, G. Gibert. Dependency of slab geometry on absolute velocities and conditions for cyclicity: insights from numerical modelling. Geophysical Journal International, 2012, 189 (2), pp.747-760. 10.1111/j.1365-246X.2012.05426.x. hal-00715971

HAL Id: hal-00715971

https://hal.science/hal-00715971

Submitted on 26 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dependency of slab geometry on absolute velocities and conditions for cyclicity: insights from numerical modelling

G. Gibert, M. Gerbault, R. Hassani and E. Tric

Université Nice-Sophia Antipolis, Centre National de la Recherche Scientifique (C.N.R.S.), Institut de Recherche pour le Développement (I.R.D.), Laboratoire GéoAzur, 250 rue Albert Einstein, 06560 Valbonne, France. E-mail: hassani@unice.fr

Accepted 2012 February 20. Received 2012 February 9; in original form 2011 October 8

## SUMMARY

The aim of this study is to quantify the relationship between the kinematics of subduction, deformation in the overriding plate and the evolution of slab geometry. A 2-D finite element numerical code is used, and a first objective consists in benchmarking previously published analogue models. Far-field plate velocities are applied, and once the subducting plate reaches the 660 km discontinuity, modelled as a rigid base, we obtain two different forms or styles of subduction that depend on the overriding plate velocity $v_{op}$: if $v_{op} > 0$, the slab lies forwards on the 660 km discontinuity (style 1), and if $v_{op} \le 0$, the slab lies backwards on the discontinuity (style 2). We also obtain a cyclic pattern with the slab folding on itself repeatedly when $v_{sp} > 0$ and $2v_{op} + v_{sp} > 0$ (where $v_{sp}$ is the subducting plate velocity). These conditions result from the analysis of several simulations in which the subduction velocities and plate viscosities are varied. When the slab periodically folds on the 660 km discontinuity, periods of shallow slab dip and compression in the overriding plate are followed by periods of slab steepening and relative extension in the overriding plate. Folding periodicity is controlled by the slab viscosity and subduction velocity. When a low-viscosity zone is included in the overriding plate, the trench motion is effectively decoupled from the overriding plate velocity, therefore allowing it to be directly controlled by the deep dynamics of the slab. For the cyclic style 2 corresponding to forward folding of the slab, the low-viscosity region in the overriding plate increases the stress amplitudes oscillations, the trench motion and the folding periodicity with time. Therefore the strength of the entire overriding plate is shown to directly control the dynamics of subduction. Using the Nazca and South American plate velocities we produce models of cyclic folding with a period of ca. 22 Ma and a minimal dip angle of ca. 10°. Episodic folding of the slab on the 660 km discontinuity would produce the necessary changes in slab dip and overriding plate deformation that explain episodes of volcanic quiescence alternating with greater rates of shortening along the Andes.

**Key words:** Subduction zone processes; Dynamics of lithosphere and mantle; Kinematics of crustal and mantle deformation; Mechanics, theory, and modelling.

## 1 INTRODUCTION

Tomographic images show that the behaviour of slabs that subduct into the mantle is highly variable, with a significant number of slabs (Japan, Aleutian and South American subduction zones) that stagnate in the transition zone around 660 km depth (e.g. Fukao *et al.* 2001), while other subducting plates appear to descend deeply into the lower mantle (Fukao *et al.* 2009). These observed differences suggest that the interaction of the slab and the 660 km discontinuity influences the dynamics of subduction zones, the temporal evolution of upper plate deformation and slab dip. On the other hand, migration of volcanic arcs over time is often related to changes in slab dip. In the Southern and Central Andes, for instance, three periods of 30–40 Ma of east–west migration of the magmatic arc as-

sociated to variations in continental shortening have been attributed to changing slab dip (Haschke *et al.* 2002; Ramos 2009). While this cyclic property of Andean magmatism and orogenic growth has been interpreted as resulting from periodic continental delamination (Kay & Mpodozis 2002; DeCelles *et al.* 2009; Ramos 2009), the influence of deeper processes occurring in the mantle is also a possibility.

Numerous modelling studies have been devoted to the dynamics of subduction zones, to quantify the links between plate tectonics on the Earth's surface (plate motions, trench migration, backarc regimes of deformation), the slab geometries and the coupling with mantle flow. The approach used is either statistical, numerical or analogue. Each one of these has advantages and disadvantages, although they are all complementary. Nevertheless, studies of the
temporal evolution of the slab and its influence on upper plate de-
formation can only be carried out starting from the numerical or
the analogue approaches. Their principal drawback is linked to the
choice of the boundary conditions. Assumptions about the appropri-
ate boundary conditions for subduction dynamics have been debated
for several decades, and two classes of approaches have emerged
over the years. The first class of approaches adopts the point of
view of mantle convection and considers the predominant role of
the Stokes velocity of descending slabs into the flowing mantle. The
opposing approach is based on the point of view of plate tecton-
ics, and accounts for the kinematic constraints imposed at the top
surface of the Earth.

Both modelling viewpoints have already yielded important an-
swers to subduction dynamics. Numerical models based on as-
sumptions from mantle dynamics have explored the interaction of
a freely descending slab into the mantle. The viscosity and den-
sity ratios of a slab with the surrounding mantle were identified as
the key parameters, and determine its Stokes velocity (Christensen
1996; Ciskova et al. 2002; Funiciello et al. 2003b; Schellart 2005;
Goes et al. 2008; Ribe 2010). The shape of a slab as it falls upon
the 660 km depth discontinuity depends on this Stokes velocity
and its relative stiffness, producing geometries of either backward
draping with a recumbant fold, or folding and piling up forwards
(Royden & Husson 2006; Van Hunen et al. 2008; Ribe 2010;
Schellart 2010). Schellart (2005) noted that during these draping
and folding processes, episodic hinge-migration of the slab oc-
curred at the surface. Stegman et al. (2010a,b) and Capitanio et al.
(2010) produced numerically a variety of 3-D slab geometries, or
styles, some of which correspond to those modelled in the labora-
tory by Schellart (2005), Heuret et al. (2007) and Guillaume et al.
(2009). In these approaches trench motion was used as a key fac-
tor that relates the styles of subduction and deformation at the top
surface.

In relation to slab anchoring conditions at depth, steep slabs are
understood to induce minor shortening or extension in the overriding
plate, as opposed to shallowly dipping slabs that induce greater
shortening of the overriding plate (Jarrard 1986; Ciskova et al. 2002;
Funiciello et al. 2004; Schellart 2005). Conversely, statistical data
and analogue modelling show that the overriding plate velocity con-
trols two extreme styles of slab geometry, lying either forwards or
backwards on the 660 km depth discontinuity (Heuret & Lallemand
2005). These results indicate that plate kinematics imposed by spa-
tial constraints on the Earth's surface have a role in slab dynamics,
which brings us to the second modelling point of view. In this sec-
ond approach, rather using assumptions from solid-mechanics, the
kinematics and the strength of tectonic plates, together with inter-
plate friction were shown to have first-order effects on the arc and
backarc stress regimes (Hassani et al. 1997; Chemenda et al. 2000;
Gerbault et al. 2009; Yamato et al. 2009).

When using plate kinematics as input for boundary conditions,
non-stationary behaviour of the slab associated with tectonic pulses
within the overriding plate was an important achievement that was
obtained with analogue models by Guillaume et al. (2009). How-
ever, the validation of this process has not been clearly demonstrated
numerically yet because it remains technically delicate to account
altogether for the free-surface, the elasto-plastic solid behaviour
of tectonic plates, and the viscous flow in the upper mantle, in a
fully self-consistent way. Numerically, cyclic stress regimes in the
overriding plate were obtained when free subduction (e.g. weight-
driven) and a kinematically fixed overriding plate were considered
(Capitanio et al. 2010). However, as will be discussed later on, these
models neglect the potential strength of the plates, and do not con-

sider the Earth's surface stress field and kinematics as an integral
part of the dynamics of subduction.

In this study, we will adopt the second modelling point of view.
In Section 2, the mechanical problem and the numerical approach
are presented in detail. In Section 3, we benchmark the results from
the analogue modelling (Heuret et al. 2007; Guillaume et al. 2009),
and reproduce the evolution of the dip of a slab as it falls over a rigid
660 km discontinuity. Once this benchmark is achieved, we aim at
establishing a relationship between the kinematics of a subduction
zone, the slab geometry and the regime of deformation of the over-
riding plate, in comparison to Funiciello et al. (2003a), Guillaume
et al. (2009) and Heuret et al. (2007). A kinematic phase diagram is
presented in Section 4, that demonstrates the conditions for cyclic
styles of subduction, deduced from a series of numerical tests that
show the impact of absolute plate velocity on the morphology of
the slab down to the 660 km discontinuity. We propose a velocity
criterion from which the subduction zone evolves with time, which
produces either steady-state or cyclic styles and controls both the
shortening regime of the overriding plate and the slab geometry.
We also show how the geometry of the subduction is sensitive to
the slab viscosity, which is in agreement with previous studies, but
more importantly, how the trench motion is sensitive to the strength
of the overriding plate. In Section 5 we analyse and discuss our
results with respect to other studies and data.

2 MECHANICAL MODELLING

## 2.1 Model set-up and general assumptions

We use the same approach as that presented in Hassani *et al.* (1997) and Bonnardot *et al.* (2008a,b), which is similar to that used by Buiter *et al.* (2001). The model consists of two initially horizontal lithospheric plates of density ρ₁ separated by a pre-defined fault and overlying an asthenospheric mantle of density ρₘ (see Fig. 1). To reproduce the analogue modelling results presented in Guillaume *et al.* (2009) and Heuret *et al.* (2007) we first run numerical experiments under laboratory conditions and scales. Temperature effects are not taken into account and terrestrial curvature is neglected. As in these analogue models, the interface between the upper mantle and the lower mantle is viewed as an impassable barrier where the slab remains anchored (no slip allowed). Moreover, as in Bonnardot *et al.* (2008b) and Hassani *et al.* (1997), the viscosity of the surrounding upper mantle is neglected, assuming a high-viscosity contrast between the lithosphere and upper mantle. This strong assumption was found to be consistent with the viscosity contrast used in the analogue experiments where the viscosity ratio between the subducting plate and the upper mantle is fixed to 12 000 in Heuret *et al.* (2007) and varies between 6000 and 15 000 in the study of Guillaume *et al.* (2009). We will discuss later the lower viscosity contrasts as considered by other authors (Schellart 2009; Capitanio *et al.* 2010). At first, velocity boundary conditions in the numerical model are set equal to those in the analogue model, and then they are varied to test their effects on the subduction dynamics.

The approach that we take (the second viewpoint) has previously
been referred to as 'kinematic' models of subduction because plate
motions are imposed; however, they are still dynamic models, as the
slab deformation and stress distribution are self-consistent results of
the model. The other approaches based on the slab Stokes velocity
(e.g. Capitanio et al. 2010; Stegman et al. 2010a) also imply some
kinematical artefacts because: (i) the top surface is not free; and/or
(ii) the overriding plate, if present, is either fixed or set free in the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic representation of the model at the initial time with the geometrical and mechanical parameters and boundary conditions. The lithospheric plates are defined as a viscoelastic Maxwell body of viscosity η and density ρ_l, while the upper mantle is an non-viscid passive fluid of density ρ_m. The upper mantle-lower mantle discontinuity is modelled as an impermeable barrier on which the slab will be anchored due to a sufficiently high friction coefficient. Note that velocities are expressed in the framework of this discontinuity and that positive values of the normal velocities component, v_op and v_sp, correspond to trenchward motion.</figcaption></figure>

horizontal direction. Our ‘kinematic’ models differ from such ‘free-subduction’ models because we account for the so-called ‘external forces’ on the dynamics of subduction (Goes et al. 2008), which we regard instead as an integral part of these dynamics. This view is validated by the observation by Lallemand et al. (2008) of a poor statistical correlation between the overriding plate velocity and the subducting plate velocity, with a quality factor R² < 0.37. These authors also commented that the overriding plate velocity is twice as influential on the tectonic regime of the upper plate as is the subducting plate velocity. In this study, we will make the point here that the spatial constraints imposed by the Earth’s surface indeed have a first-order role in the dynamics of subduction.

## 2.2 Governing equations

The governing equations of the quasi-static evolution of lithospheric plates occupying at time t the physical domain Ω_t ⊂ ℝ² are given by

$$ \left\{ \begin{array}{l} \mathrm{div}\boldsymbol{\sigma} + \rho_l \boldsymbol{g} = \mathbf{0} \quad \text{in } \Omega_t, \\ \frac{D\boldsymbol{\sigma}}{Dt} = \mathcal{M}(\boldsymbol{\sigma}, \mathbf{d}) \quad \text{in } \Omega_t, \end{array} \right. \qquad (1) $$

where σ is the Cauchy stress field, g is the vector of gravity acceleration, ρ_l is the lithospheric density and d = 1/2(∇v + ∇v^T) is the symmetric gradient of the velocity field v. The symbol D/Dt stands for an objective time derivative (Jaumman or Green-Naghdi derivative) and the function M represents a general constitutive law. To make the comparison between analogue and numerical results easier, the M-laws used in this study correspond to the elastic and the Maxwell viscoelastic rheologies:

$$ \mathcal{M}(\boldsymbol{\sigma}, \mathbf{d}) = 2G\mathbf{d} + \lambda\mathrm{tr}(\mathbf{d})\mathbf{I} + \frac{G}{\eta}\mathrm{dev}\boldsymbol{\sigma}, \qquad (2) $$

where I is the identity tensor, tr and dev are the trace and deviatoric operators, respectively, λ and G are the Lamé parameters and η is the viscosity (η = ∞ for a purely elastic material).

Unilateral constraints are taken into account on the contact surface Γ (subduction interface) between the two plates. These

constraints are given by the Signorini relation (condition of no interpenetration) and the Coulomb friction law

$$ \left\{ \begin{array}{l} \delta v_n \le 0, \quad \sigma_n \le 0, \quad \delta v_n \sigma_n = 0, \\ |\sigma_t| \le -\mu \sigma_n \quad \text{if } \delta v_t = 0, \\ \sigma_t = \mu \sigma_n \frac{\delta v_t}{|\delta v_t|} \quad \text{if } \delta v_t \ne 0, \end{array} \right. \qquad (3) $$

where δv_n(x) and δv_t(x) are, respectively, the normal and tangential components of the relative velocity between the plates at a point x, μ is the effective Coulomb friction coefficient and σ_n(x) and σ_t(x) are the normal and tangential stresses. As in analogue modelling, the 660 km discontinuity is simulated as an infinitely rigid plane. Thus, identical frictional conditions (3) rule the interaction between the slab and this discontinuity.

## 2.3 Initial geometry and boundary conditions

The plate length is set to 3000 km to make the effects of the border negligible. We assume no initial topography on the overriding plate. The two plates are simply pushed against each other by far-field boundary conditions (indices ‘op’ and ‘sp’ stand for ‘overriding plate’ and ‘subducting plate’, respectively):

$$ \mathbf{v} \cdot \mathbf{n} = -v_{\mathrm{op}} \quad \text{on } \Gamma_{\mathrm{op}} \quad \text{and} \quad \mathbf{v} \cdot \mathbf{n} = -v_{\mathrm{sp}} \quad \text{on } \Gamma_{\mathrm{sp}}, \qquad (4) $$

where v is the velocity vector and n is the outward normal vector of the vertical edges Γ_op and Γ_sp of the two plates (Fig. 1). Because of the negative sign in (4), positive values of v_op and v_sp correspond to trenchward motion.

Based on the analogue experiments of Heuret et al. (2007) and Guillaume et al. (2009), in which a high viscosity ratio between the plates and the upper mantle is used, and on the numerical study of Bonnardot et al. (2008a), where it was shown that viscous effects can be neglected if the upper-mantle viscosity is lower than or equal to 10¹⁹ Pa s, we infer that the plates/upper-mantle interaction can be reduced to its isotropic part. Each interface, Γ_m, of the lithospheric plate in contact with the asthenosphere is then subjected to a hydrostatic pressure:

$$ \mathbf{n} \cdot \boldsymbol{\sigma} \mathbf{n} = P_m := \frac{1}{\beta} \ln(1 - \beta \rho_m^0 gz) \quad \text{on } \Gamma_m, \qquad (5) $$
where $\beta$ is the compressibility modulus of the upper mantle (assumed constant), $\rho_m^0$ is its density at the base of the lithosphere, $z$ is the distance to the hydrostatic level and $g = \|g\|$ (see Hassani et al. 1997; Bonnardot et al. 2008b, for more details).

An approximate solution to the continuous problem (1–5) is computed using the finite element code ADELI. ADELI has been used in numerous geodynamical applications, for processes at the crustal scale (e.g. Vanbrabant et al. 1997; Huc et al. 1998; Berger et al. 2004; Got et al. 2008) as well as at the lithospheric scale (e.g. Hassani et al. 1997; Lesne et al. 2000; Chéry et al. 2001; Bonnardot et al. 2008a,b; Neves et al. 2008). This code belongs to the FLAC family of codes (Cundall & Board 1988; Poliakov & Podladchikov 1992), and is based on the dynamic relaxation method (Underwood 1983). Numerical aspects about ADELI in two and in three dimensions can be found in Hassani et al. (1997) and Chéry et al. (2001).

## 2.4 Description of the reference analogue modelling

Due to the difficulty of determining analytical solutions, the best way to validate a dynamic process of subduction is to benchmark laboratory experiments with a numerical approach. To this end, we briefly summarize the experimental results of Heuret et al. (2007) and Guillaume et al. (2009). Then, one of these analogue experiments is reproduced numerically and is followed by a parametric study that tests numerically the impact of plate absolute velocities and viscosities upon the slab geometry and the overriding plate deformation.

Heuret et al. (2007) studied a large set of models with various relative subduction velocities. For each relative plate motion, several pairs of overriding and subducting plate velocities were used. In their experiments, Guillaume et al. (2009) used different plate thicknesses, viscosities, densities and subduction velocities, while keeping the overriding plate fixed. Both analogue experiment set-ups were composed of viscoelastic silicone plates (lithosphere) placed above glucose syrup (mantle) in a Plexiglas tank. The bottom of the Plexiglas tank represented the 660 km discontinuity

between the upper mantle and the lower mantle. Velocities were imposed by a rigid piston moving towards or away from the trench. With their simulations, Heuret et al. (2007) identified two styles of subduction, depending on the overriding plate velocity direction: when the overriding plate was moving towards the trench, the slab folded along the bottom of the box [named style 1 after Heuret et al. (2007), fig. 4a of their article] and the overriding plate shortened. When the overriding plate was moving away from the trench, the slab bent backwards (named style 2, fig. 4b of Heuret et al. 2007), and the overriding plate remained nearly undeformed.

Guillaume et al. (2009) observed the stacking of the slab on the bottom of their modelling box (due to the fixed trench position and the anchoring of the slab). This stacking of the slab results in periods of flattening and steepening of the dip angle, which correspond to episodes of increased and decreased shortening in the upper plate. These experimental results show that different regimes (cycles or no cycles) can occur while the kinematic boundary conditions remain constant. This result has strong implications when reconstructing plate motion from field observations of shortened tectonics units. Our first goal here is to reproduce these results and then to identify the mechanical conditions that favour the occurrence of cyclic slab dip.

# 3 RESULTS

## 3.1 Comparison of numerical and analogue experiments

We first use the same reference model as in Guillaume et al. (2009) and adopt the laboratory scale (see Table 1 for the model parameters). We check that two numerical experiments give exactly the same results as long as the values of the following dimensionless numbers remain the same:

$$ 
\begin{aligned}
\kappa_1 &= \frac{\eta v_s}{E e}, & \kappa_2 &= \frac{\Delta \rho g e}{E}, & t^* &= \frac{\eta}{E T}, & h^* &= \frac{h}{e}, \\
v_{op}^* &= \frac{v_{op}}{v_s}, & v_{sp}^* &= \frac{v_{sp}}{v_s}, & &
\end{aligned}
\quad (6) $$

<table><caption>**Table 1.** Mechanical parameters used for the comparison test between the numerical modelling and the analogue models of Guillaume et al. (2009). The underlined values are those used by Guillaume et al. (2009) in their reference model. The other parameters ($E$, $v$, $\mu$) are not well known in the analogue set-up and their values were tested numerically. Note that the viscosity ratio $\eta_m/\eta$ between the upper mantle and the plate is non-zero in the analogue experiments, but very small ($1.6 \times 10^{-4}$). Corresponding values at the natural scale are indicated for reference only, in the third column. These values are computed using the similarity criteria given by (6) and by prescribing that $\Delta\rho_{\text{nature}} = \Delta\rho_{\text{model}}$, $v_{\text{nature}} = 10 \text{ cm a}^{-1}$ and $e_{\text{nature}} = 90 \text{ km}$.</caption><thead><tr><th>Parameters</th><th>Values used in the laboratory scale model</th><th>Corresponding values at the natural scale</th></tr></thead><tbody><tr><td>Plate thickness, <i>e</i></td><td><u>13 mm</u></td><td>90 km</td></tr><tr><td>density contrast, Δρ = ρ<sub>1</sub> - ρ<sub>m</sub></td><td><u>76 kg m<sup>-3</sup></u></td><td>76 kg m<sup>-3</sup></td></tr><tr><td>Subducting plate velocity, <i>v</i><sub>sp</sub></td><td><u>1 cm min<sup>-1</sup></u></td><td>10 cm a<sup>-1</sup></td></tr><tr><td>Simulation time, <i>T</i></td><td><u>3300 s</u></td><td>38.8 Ma</td></tr><tr><td>Plate viscosity, η</td><td></td><td></td></tr><tr><td>Subducting plate</td><td><u>5 × 10<sup>5</sup> Pa s</u></td><td>1.28 × 10<sup>24</sup> Pa s</td></tr><tr><td>Overriding plate</td><td><u>3 × 10<sup>5</sup> Pa s</u></td><td>7.7 × 10<sup>23</sup> Pa s</td></tr><tr><td>Young modulus, <i>E</i></td><td>7000 Pa</td><td>4.85 × 10<sup>10</sup> Pa</td></tr><tr><td>Poisson ratio, <i>ν</i></td><td>0.25</td><td>0.25</td></tr><tr><td>Friction coefficients, μ</td><td></td><td></td></tr><tr><td>Plate/plate</td><td>0.015</td><td>0.015</td></tr><tr><td>Plate/discontinuity</td><td>0.095</td><td>0.095</td></tr><tr><td>Viscosity ratio, η<sub>m</sub>/η</td><td>0</td><td>0</td></tr></tbody></table>
where $v_s = v_{sp} + v_{op}$ is the subduction velocity ($v_s = v_{sp}$ in the reference model), $T$ is the duration of the experiment, $e$ is the plate thickness, $E$ is the Young modulus, $\Delta\rho = \rho_l - \rho_m$ is the density contrast between the lithosphere and upper mantle, $\eta$ is the plate viscosity and $h$ is the model thickness (see Fig. 1). The corresponding parameters at the natural scale are given in Table 1. Note that three physical properties are not well constrained in the analogue modelling: the Young's modulus of the silicone putty (used to model the lithosphere), the equivalent friction coefficient of the silicone/silicone contact and that of the silicone/box (plexiglas), and the initial fault dip. A Young's modulus of 5–7 kPa (which is, to our knowledge, an acceptable range for silicone putty) and a 30° fault-dip angle were assumed in our numerical experiment. Friction between the bottom of the box and the slab was set to a value that

only allows minor displacement of the slab, as is the case in the analogue models.

With this chosen set of values and despite the inviscid nature of our upper mantle, the results at different time steps show slab geometries remarkably close to those observed in the analogue models (Fig. 2). Note however the larger amount of shortening that the overriding plate undergoes in the analogue model (which is particularly obvious at the end of the experiment). This disagreement can be explained either by a difference in rheology, which makes the overriding plate weaker in the experimental set-up than in the numerical model, or by the use of the Coulomb contact friction law, which might not be well adapted to model numerically the silicone/silicone contact and silicone/plexiglas contact, especially in the presence of a liquid. Indeed, in the analogue experiment, the

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison between numerical result (a) obtained for the reference case, and the analogue result of Guillaume *et al.* (2009) (b). The general evolution is very similar for the two approaches. There are however two major differences: (1) the trench position is almost unchanged in (a) while advancing in (b); (2) the deep slab is anchored in (a), while it slips on the box bottom at the end of the experiment in (b). Frictional coupling between the two plates and between the slab and the box bottom, which is modelled by using constant Coulomb friction coefficients in (a), probably evolves during the experiment in (b).</figcaption></figure>
<table><caption>**Table 2.** Mechanical parameters used in our numerical experiments.</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Plate thickness, <i>e</i></td><td>90 km</td></tr><tr><td>Density contrast, Δρ = ρ<sub>1</sub> – ρ<sub>m</sub></td><td>50 kg m<sup>-3</sup></td></tr><tr><td>Subducting and overriding plate viscosities, η</td><td>10<sup>22</sup> to 10<sup>26</sup> Pa s</td></tr><tr><td>Young modulus, <i>E</i></td><td>5 × 10<sup>10</sup> Pa</td></tr><tr><td>Poisson ratio, <i>ν</i></td><td>0.25</td></tr><tr><td>Friction coefficients, μ</td><td></td></tr><tr><td>Plate/plate</td><td>0.01</td></tr><tr><td>Plate/discontinuity</td><td>0.2</td></tr><tr><td>Viscosity ratio, η<sub>m</sub>/η</td><td>0</td></tr><tr><td>Subducting and overriding plate velocities, <i>v</i><sub>op</sub>, <i>v</i><sub>sp</sub></td><td>–20 to 20 cm a<sup>-1</sup></td></tr></tbody></table>

overriding and subducting plates are initially decoupled by a thin layer of glucose syrup, but the thickness of this layer decreases during the convergence, thus the interplate coupling should increase with time and consequently, intensify the overriding plate shortening. We also note that at the end of the experiment, the slab is perfectly anchored at the base in the numerical model, whereas it undergoes relatively important slip at the base of the analogue model box. In summary and to a first order, our numerical simulations are comparable to analogue models.

## 3.2 Effects of absolute plate velocities

Since subduction is controlled by far-field boundary conditions in our models, and because the overriding plate deformation remains small (except for low-viscosity upper plates; see Section 3.4), the absolute trench velocity $v_t$ is almost equal to the overriding plate velocity $v_{op}$. Thus the subduction velocity (relative convergence velocity) is simply given by $v_s = v_{op} + v_{sp}$. As already mentioned, positive values of $v_{op}$ or $v_{sp}$ represent motion towards the trench (see Fig. 1), and obviously, only positive values of $v_s$ are considered (so that subduction takes place).

To understand the impact of absolute velocities on the slab geometry, we performed four sets of tests that considered four fixed values of $v_s$ (16 experiments). Thus, each test corresponded to a different pair ($v_{sp}$, $v_{op}$) (a sampling point of the half-plane $v_s > 0$), and all of the other parameters were held constant with tectonic-scale values (see Table 2).

The resulting slab geometries are represented in a ($v_{sp}$, $v_{op}$) diagram in Fig. 3(a) where we make use of the same graphical organization used by Heuret et al. (2007). On this graph, results that correspond to each of the four sets of experiments at constant relative convergence velocity belong to the same diagonal line.

Every experiment starts in a similar manner, and the slab descends through the asthenosphere with a dip that depends on the subduction

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Relationship between absolute velocities $v_{op}$ and $v_{sp}$ and slab geometry: (a) Figure organization is based on Heuret et al. (2007). Each cross-section corresponds to a pair ($v_{sp}$, $v_{op}$) represented by dots; each diagonal line corresponds to results obtained with constant $v_s$. To facilitate the comparison with Heuret et al. (2007), the same velocity unit (cm min$^{-1}$) is used in this diagram. Thus, 1 cm min$^{-1}$ corresponds to 10 cm a$^{-1}$ at the natural scale. The red zone, bounded by the lines $v_{sp} = 0$ and $2v_{op} + v_{sp} = 0$, represents the pairs of absolute velocity for which the slabs show cyclic behaviour. (b) The plate velocities are normalized by the subduction velocity $v_s = v_{op} + v_{sp}$. Cycles arise when $|v_{op}^*| < 1$, while the stretching mode occurs when $v_{op}^* > 1$ for style 1, and $v_{op}^* < -1$ for style 2.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Four characteristic modes of subduction. The colour code corresponds to the second invariant of deviatoric strain, and the arrows show far-field absolute plates velocities (black arrows indicate trenchward motion). The symbol ⊗ indicates that the velocity is zero.</figcaption></figure>

velocity ($v_s$). However, once it reaches the 660 km transition zone and anchors to it, the way in which it deforms depends strongly on the individual values of plate velocities ($v_{sp}$ and $v_{op}$).

As shown in Fig. 3(a), four distinguishable modes of subduction occur, and these are shown in Fig. 4: the slab can lie forwards (style 1) on the 660 km discontinuity, where it can fold episodically (mode 1a), or it can stretch (mode 1b). Alternatively, the slab can lie backwards (style 2) and cycles of folding can take place (mode 2a) or not (mode 2b).

The subduction modes are organized logically as seen in Fig. 3(a). First, conditions for the development of styles 1 and 2 are obvious: the slab lies forwards for positive upper plate velocities ($v_{op} > 0$), that is, when the upper plate moves towards the trench. In contrast, the slab lies backwards when the upper plate is fixed or retreating ($v_{op} \le 0$).

These conditions have already been established from the analogue results of Heuret et al. (2007) and Guillaume et al. (2009). However, the cyclic regime for mode 1a was not observed by either Heuret et al. (2007) or Guillaume et al. (2009), and the only case studied by Guillaume et al. (2009) that displayed periodic slab folding corresponds to $v_{op} = 0$.

As shown in Fig. 3(a), cyclic behaviour takes place only when $v_s > |v_{op}|$, that is, when the absolute subduction velocity is greater than the upper plate velocity. Note that this condition is equivalently expressed by $v_{sp} > 0$ and $2v_{op} + v_{sp} > 0$ ($v_s = v_{op} + v_{sp}$ is assumed positive for subduction to occur).

The conditions for which each of the four modes hold can therefore be specified as follows:

Style 1: the slab lies forwards if $v_{op} > 0$ and

(i) mode 1a: the slab undergoes folding when $v_{sp} > 0$,

(ii) mode 1b: the slab can stretch (elongates) when $v_{sp} \le 0$,

Style 2: the slab lies backwards if $v_{op} \le 0$ and

(i) mode 2a: the slab folds when $v_{sp} + 2v_{op} > 0$,

(ii) mode 2b: the slab can stretch when $v_{sp} + 2v_{op} \le 0$.

We also note that perfect steady-state regimes (neither folding nor stretching of the slab) are only those for which $v_s = |v_{op}|$, that is $v_{sp} = 0$ for style 1 or $2v_{op} + v_{sp} = 0$ for style 2. Thus, in

Fig. 3(a) steady-state regimes correspond to the boundary of the red zone that delimits the occurrence of subduction with and without cycles.

Another way to describe this result is to make use of the normalized velocities $v_{op}^* = v_{op}/v_s$ and $v_{sp}^* = v_{sp}/v_s$, from which it is straightforward to determine that the subduction styles fall on the line $v_{sp}^* + v_{op}^* = 1$ and that (see Fig. 3b):

(i) a cyclic regime occurs when $-1 < v_{op}^* < 1$ (mode 1a if $v_{op}^* > 0$, otherwise mode 2a),

(ii) when $v_{op}^* > 1$ the slab stretches according to mode 1b,

(iii) when $v_{op}^* < -1$ the slab stretches according to mode 2b,

(iv) steady-state regimes occur for $|v_{op}^*| = 1$ (style 1 if $v_{op}^* = 1$, style 2 if $v_{op}^* = -1$).

Depending on the slab viscosity, the occurrence of folds increases (folds are closer to each other) when the ratio $v_{op}/v_s$ decreases, that is, when the contribution of the overriding plate velocity to the subduction velocity decreases, as shown in Fig. 5.

### 3.3 Effects of slab viscosity

To illustrate the role of slab viscosity, $\eta$, and to compare the results with other studies (see the discussion in Section 4.3), we performed another set of tests by varying $\eta$ between $10^{22}$ Pa s and $10^{26}$ Pa s. We consider that the value of $10^{24}$ Pa s is a most realistic average, as has also been considered in many studies (Heuret & Lallemand 2005; Royden & Husson 2006; Van Hunen et al. 2008; Guillaume et al. 2009). Each value of the viscosity is investigated with four velocity pairs that correspond to each previously highlighted subduction modes.

When $\eta = 10^{26}$ Pa s, the model corresponds to the extreme case of a viscoelastic slab behaving purely elastically for the time considered. This high viscosity yields a high stiffness, and thus less bending associated to the descent of the slab into the mantle (Ribe 2010). The less viscous cases ($\eta = 10^{22}-10^{23}$ Pa s) produce slab dripping as described in previous studies (e.g. Christensen 1996; Ribe 2010). At the other extreme, low slab viscosities can even lead to slab detachment, depending on the value of $v_{sp}$ versus the slab's Stokes velocity (Billen 2008; Ribe 2010; Duretz et al. 2011).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The period of the cycles in mode 1a depends directly on the normalized velocity $v_{op}^*$. The solution corresponds to four values of $v_{op}^*$ and with $\eta = 10^{24}$ Pa s.</figcaption></figure>

When $\eta = 10^{24}$ Pa s as in the reference model (Section 3.1), mode 1a cyclic folding adopts a variable slab dip with depth and time. When the slab forms a fold, it adopts a low dip between the base of the lithosphere and ca. 300 km depth, and a dip >90° at greater depths. When the fold lies down, the slab recovers the dip it had prior to reaching the discontinuity, and a new cycle starts (Fig. 6). When varying slab viscosity, this process is affected in the main following ways (Fig. 6):

(1) The dip of the slab as it reaches the 660 km discontinuity decreases with increasing viscosity, because resistance to bending is greater. Folds can still develop, but are limited by the maximum vertical size of the domain (Royden & Husson 2006; Ribe 2010). Capitanio et al. (2010) did not obtain such a significant effect of slab viscosity on the geometry of subduction, simply because their viscosities were too small (see discussion in 4.3). As shown by Ribe (2010) and Royden & Husson (2006), slab stiffness starts to have a significant effect for viscosities greater than $10^{24}$ Pa s.

(2) Slab viscosity does not modify the subduction style, in the sense that conditions determined for the orientation of deposition and the formation of folds remain the same (Fig. 6a). However, increasing viscosity (and thus slab stiffness) increases the length over which the slab bends when forming a fold over the kilometre discontinuity (Ribe 2010). Thus the fold forms, but over a longer time period (Fig. 6b). In other words, cycle duration increases with increasing viscosity, because of a greater bending moment (and greater time for the folds to develop).

### 3.4 The role of a weak zone in the overriding plate

Here, to assess the impact of the strength of the overriding plate on the dynamics of subduction and trench motion, we performed two additional tests in which a portion of the overriding plate behaves as a weaker viscoelastic domain. Initial conditions are chosen with $v_{op} = 12.6 \text{ cm a}^{-1}$ ($4 \times 10^{-9} \text{ m s}^{-1}$), and $v_{sp} = 6.3 \text{ cm a}^{-1}$ ($2 \times 10^{-9} \text{ m s}^{-1}$), which puts us within the previously defined kinematical mode 1a, with cyclic folding forwards on the 660 km depth boundary.

We visualize the results by comparing variations in the horizontal length ($L$) of the weak zone (Fig. 7). In the reference model with a uniform and stiff overriding plate (plate viscosities were set equal to $10^{24}$ Pa s), the trench motion $v_t$ almost equalized the motion of the overriding plate, $v_t \simeq v_{op}$. We observed a minor effect of the slab cyclic folding on $v_t$ (black circles in Fig. 7). Due to the slab-pull effect and the low interplate friction, the upper plate undergoes stretching, as shown by the positive slope of the linear portion of this curve. Variations around this tendency (about $\pm 1$ per cent) show however that periods of shortening alternate between periods of stretching.

If the viscoelastic block is assigned a viscosity of $5 \times 10^{23}$ Pa s, half that of the subducting plate, then this portion of the overriding plate acts as an inelastic spring under a lower state of stress. Trench motion becomes mechanically disconnected to the overriding plate, and conversely, oscillates in response to the slab behaviour at 660 km depth (blue triangles in Fig. 7).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Effects of slab viscosity. (a) The slab geometry (lower edge of the slab) for different viscosities and at different times shows that folds occur roughly at the same place [except for the low-viscosity case ($\eta = 10^{22}$ Pa s)]. (b) As for (a), but plotted only for the first fold. The triangular symbols mark the trench location. (c) Successive slab geometries from the beginning of the formation of two consecutive folds (blue lines) for $\eta = 10^{24}$ Pa s, to its end.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Effects of a weak viscoelastic zone (grey area) in the upper plate. Variations in length ($L$) of the weak zone with time, (a) for the reference case with a homogeneous viscoelastic plate ($\eta_{weak} = \eta = 10^{24}$ Pa s), (b) for a weak zone half as viscous as the remainder of the plate ($\eta_{weak} = 5 \times 10^{23}$ Pa s) and (c) for a viscoelastic weak zone five times less viscous ($\eta_{weak} = 2 \times 10^{23}$ Pa s).</figcaption></figure>

If the viscoelastic block is assigned a viscosity five times lower than that of the subducting plate ($2 \times 10^{23}$ Pa s), then oscillations in trench motion occur with even greater amplitudes (red diamonds in Fig. 7). Folding-induced oscillations can contribute to nearly half of the variations in total length of the overriding plate (15 per cent of the initial weak zone length $L$). The periodicity of these oscillations increases with time, together with the periodicity of the folds on

the 660 km discontinuity. This increase in periodicity is linked to a progressive thinning of the overriding plate, down to about half its original thickness. Consequently its capacity to store stress is reduced, which in turn affects the slab by increasing its folding periodicity over the 660 km depth boundary.

Further decrease in the overriding plate viscosity will approach the idealistic situation of no overriding plate, thus rendering trench
motion only dependent on the slab subduction velocity ($v_t \simeq v_s$). This was shown by ‘free subduction’ models (Schellart 2008; Van Hunen et al. 2008). While ‘free-subduction’ models showed that slab stiffness is the primary control on trench motion (Van Hunen et al. 2008), we show here that the strength of the overriding plate also has a first-order role in trench motion.

# 4 DISCUSSION

## 4.1 Analysis of modelled cyclicity for a stiff overriding plate

Our numerical results complete the slab geometry diagram proposed by Heuret et al. (2007), which relates overriding and subducting plate velocities, and lead us to propose a velocity criterion that determines the style and the behaviour of a subduction. We confirm that $v_{op}$ is the key parameter that determines the subduction style as suggested by Heuret et al. (2007). The sign of $v_{op}$ constrains the occurrence of Styles 1 and 2, with the slab taking one or the other direction when it begins to lie on the 660 km discontinuity. Until the end of the simulation, the slab dip remains approximatively constant when there is no folding. Since the friction at the 660 km discontinuity allows only little to no movement, the slab simply lies in the same direction as that of the motion of the overriding plate.

Moreover, we can infer here that the cyclic behaviour of a subducting slab is controlled by a relationship between $v_s$ and $v_{op}$. When folding occurs, the dip changes in cycles (see Fig. 6c).

In our reference configuration where the subduction velocity is imposed by far-field conditions and the overriding plate is homo-

geneous and stiff, $v_t$ is approximatively equal to $v_{op}$. The distance covered by the trench during a time $\Delta t$ represents the distance of free space along the fixed 660 km discontinuity and over which the subducted lithosphere can be deposited during that time. The subduction velocity $v_s$ determines the length of the subducting plate entering the asthenosphere during $\Delta t$. Two cases can thus be highlighted here:

(i) $|v_{op}| \simeq |v_t| \ge v_s$ and the free space on the discontinuity is long enough for the slab to lie on it,

(ii) $|v_{op}| \simeq |v_t| < v_s$ and the free space is shorter than the length of subducted lithosphere, so that cycles occur.

For this last case and recalling that $v_s = v_{op} + v_{sp}$, cycles occur in the diagram region defined by $v_{sp} > 0$ and $2v_{op} + v_{sp} > 0$ (red zone in Fig. 3a). Cycles of mode 1a are restricted to the square region, which correspond to $v_{sp} > 0$ and $v_{op} > 0$, while cycles of mode 2a belong to the triangular region defined by $2v_{op} + v_{sp} > 0$ and $v_{op} \le 0$. Cycles occur because the slab is anchored on the 660 km discontinuity. Models without any friction between the discontinuity and the slab result in steady-state processes regardless of imposed plate velocities.

As noted by Guillaume et al. (2009) for mode 2a, the slab dip angle oscillates between low and high values during subduction, in our case from 10° to 75° when cycles take place. The dip angle is a direct consequence of the interaction of the slab with the 660 km discontinuity, which imposes its geometry on the subducting plate. Cycles repeat periodically for both styles, as shown for style 1 in Fig. 8(a). As indicated by Heuret et al. (2007), the dip angle values are smaller for mode 1b than for mode 2b (cases without cycles).

<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Plot of the mean dip angle computed between 90 and 200 km depth with time for a model with the Andean velocities (Nuvel-1A): $v_{op} = 4.3 \text{ cm a}^{-1}$, $v_{sp} = 2.9 \text{ cm a}^{-1}$. (b) Far-field averaged horizontal stress in the upper plate ($< \sigma_{xx} > = \frac{1}{c} \int_0^e \sigma_{xx}|_{\Gamma_{op}} dy$) showing episodes of extensional and compressional regimes. If the frictional coupling between the two plates is increased, the plot is shifted down, and regimes of less compression alternate with regimes of greater compression.</figcaption></figure>
## 4.2 Stress and deformation in the overriding plate and trench motion

The state of stress of the overriding plate links directly with the dip of the slab (Fig. 8). Because we applied far-field plate velocities, deformation of the homogeneous overriding plate thus links directly with the values of $v_t$ and $v_{op}$, which determine the slab's motion. We observe:

(i) upper plate extension for modes 1b and 2b without cycles,

(ii) upper plate shortening when the angle of the slab shallows, and extension when it steepens for modes 1a and 2a with cycles.

In summary we obtain cyclic stress regimes in the overriding plate with constant boundary conditions.

The temporal evolution and the amount of strain in the overriding plate are directly and clearly associated with the interaction between the slab and the 660 km discontinuity. The slab is the intermediate element that controls the transfer of stress and strain between the 660 km depth and the surface. We illustrate here the fundamental role of the slab rheology on the dynamics of the system, which is complementary with previous studies that have shown that the slab age, stiffness and capacity to exchange heat with the mantle control the dynamics of subduction (Christensen 1996; King 2001; Funiciello et al. 2003b; Royden & Husson 2006; Billen 2008; Schellart 2008; Capitanio et al. 2010; Ribe 2010; Yoshioka & Naganoda 2010). We note that thermal effects are not taken into account in our study and that the slab viscosity is kept constant from the top to the bottom. In nature, the slab fluidity increases with temperature, and thus its capability to transmit stress from the 660 km depth to the surface will be reduced (e.g. Billen 2010). However, deep seismic events located around the upper-/lower-mantle discontinuity (see e.g. Fukao et al. 2009) are evidence that the slab is still elastic at depth.

When a zone of weakness is inserted into the overriding plate (modelled as a viscoelastic block in its left most extremity, Section 3.4), the trench motion $v_t$ is no longer coupled to that of the overriding plate, and it is allowed to follow the motion of the subducting plate (Fig. 7):

(i) $v_t$ can even exceed both $v_{op}$ and $v_{sp}$, because it is driven by the dynamics of folding;

(ii) $v_t$ oscillates together with alternating stretching and shortening of the viscoelastic domains;

(iii) With time, stretching of the overriding plate dominates, and oscillations in $v_t$ occur with greater amplitude. The periodicity of these oscillations also increases in time, as does the periodicity of the folds.

These results show how the strength of the overriding plate is directly involved in the dynamics of subduction, and acts as an inelastic buffer in between far-field kinematic conditions, local stress and deformation regimes and trench motion.

Topography is also linked to the slab dip angle: when the slab has a shallow dip angle, the upper plate shortens and the topography is positive. In turn, when the slab dip angle steepens, the upper plate undergoes extension and the topography becomes negative. Because the crust of the overriding plate is purely elastic or viscoelastic in our models, it is necessary to account for elastoplastic rheologies and erosion-sedimentation processes prior to any precise application to real cases. Such factors constitute a next step to future studies and are beyond the scope of this study. Nonetheless, our models might give insight into the behaviour of natural subduction zones, as is discussed in Section 4.4.

## 4.3 Slab dynamics compared to other studies

An important result of this study is that cycles are observable because $v_{op}$ and $v_{sp}$ are imposed by far-field boundary conditions on the Earth's surface, where all tectonic plates are in contact with each other and are constrained in the horizontal plane. The trench velocity is controlled by far-field boundary conditions ($v_t \simeq v_{op}$ when our overriding plate is homogeneous), in agreement with the results from Guillaume et al. (2009), Schellart (2005) and Bellahsen et al. (2005). When we incorporate a weak block in the upper plate, the trench is free to move and is controlled by the slab's motion as it deposits on the 660 km discontinuity in the easiest possible way, that is to say in a process that minimizes the number of folds.

The other important force involved in subduction dynamics is the slab pull force, which determines the Stokes velocity of the slab together with its viscosity contrast with the surrounding mantle (Schellart 2005; Capitanio et al. 2007, 2010; Goes et al. 2008; Stegman et al. 2010a; Ribe 2010). As mentioned by Goes et al. (2008), the motion induced by the slab pull (its Stokes velocity) competes with the motion imposed by far-field tectonic boundary conditions, generating either compression in the slab as it collapses and creates folds over the 660 km discontinuity, or extension with slab stretching at intermediate depths. Goes et al. (2011) attributed the departure of observed slab velocities worldwide from their modelled Stokes velocities (with freely subducting plates) to the fact that the outlawed slabs actually penetrate into the lower mantle or are controlled by basal drag or other inner-mantle processes. While these processes are important, we have made the point here that the kinematical constraints imposed by the Earth's surface are also of first-order importance, more than just exerting second-order external forces as suggested by Goes et al. (2008).

We think that there is a problem in free-subduction modelling approaches assuming homogeneous slab viscosities of about only two orders of magnitude greater than those of the surrounding mantle, with values generally much lower than $2 \times 10^{23}$ Pa s (e.g. Capitanio et al. 2007, 2010; Goes et al. 2008; Schellart 2009). A low-viscosity slab has little capacity in bending and stress storage, as demonstrated by Ribe (2010) and Royden & Husson (2006). Whereas low viscosities of the order of $10^{21}-10^{22}$ Pa s can be valid at mantle depths with increasing temperature, at the top surface where solid plate tectonics occur it remains important to consider higher viscosities. It is well known that the viscosity in the lithosphere is dominated by the power-law creep properties of olivine (Hirth & Kohlstedt 1996, 2003). For standard strain rates and geothermal gradients, the lithospheric viscosity easily exceeds $10^{25}$ Pa s above 50 km depth, with a maximum of $10^{28}$ Pa s (e.g. Arcay et al. 2006). The validity of assuming a constant plate viscosity that averages variations of nearly 10 orders of magnitude covering the top surface down to 660 km, is obviously questionable. However, this issue is beyond the scope of this study, and indeed, here too, we have assumed constant slab viscosity from top to bottom.

Energy dissipation in our models is all constrained to occur within the slab because of our inviscid mantle, as opposed to other modelling approaches in which it also occurs in the mantle (e.g. Capitanio et al. 2007). Thus in our models stress magnitudes in the slab are expected to be greater than in these other approaches, and stress transfer between 660 km depth and the surface may be slightly overestimated. For the same reason, in our models, only the slab density and viscosity act on its dip angle as it descends through the asthenosphere. When slab viscosity increases, its stiffness increases, and the slab reaches the 660 km depth discontinuity with a smaller dip (e.g. Funiciello et al. 2003b; Van Hunen et al. 2008;
Ribe 2010; Billen 2010). However, the kinematic conditions that determine the four styles of subduction that we have modelled here are not affected by such dip variations.

In contrast to our study, Capitanio *et al.* (2010) accounted for slab–mantle interactions, but they also assumed extremely weak overriding plates (a maximum strength of 35 or 58 MPa and a maximum viscosity of $3 \times 10^{21}-10^{22}$ Pa s). Such weak overriding plates can only produce small stresses, whatever the choice of boundary conditions (fixed or free upper plate). While experimental constraints on olivine for temperatures of 900°C predict yield strengths of up to 1 GPa, many studies have shown that the state of stress within tectonic plates averages to several hundred megapascals (e.g. Ord & Hobbs 1989; Hassani *et al.* 1997; Gerbault 2000; Billen 2010). We believe it is crucial to account for such strength at least in the upper part of the upper mantle. Even though, Capitanio *et al.* (2010) recognized that either fixed or free upper plate motion can affect the subduction rate ($v_s$) by up to 50 per cent (double velocities) in comparison to the other effects of plate–mantle density and viscosity contrasts. Capitanio *et al.* (2010) obtained cycles of compression and extension for their strongest upper plate when it was held fixed. We have shown here, that broader conditions for such periodicity exist, depending on both the motion and the strength of the overriding plate.

Our study is only 2-D and the third direction certainly has a strong effect on the dynamics of the system. The results obtained by Stegman *et al.* (2010b) confirm this by showing the importance of boundary conditions and slab width on trench migration velocities. Toroidal flow around slab edges acts to reduce the tractions onto the upper plate arising from corner flow circulation. In turn the incorporation of tectonic plate velocities in 3-D models are expected to yield similar first-order effects to those modeled in two dimensions here. Some of the conditions for cyclic modes of subduction that we obtain here were obtained in 3-D numerical models of free subduction, with no overriding plate, weak subducting plates, but accounting for slab–mantle interactions (Stegman *et al.* 2010a). In such models by Stegman *et al.* (2010a), the trench advances or remains quasi-stationary depending on whether the slab lies forwards or backwards. Folding modes 1a and 2a occur, depending on whether the far-field boundary of the subducting plate ($v_{sp}$) is set free or pushed towards the trench (Stegman *et al.* 2010a). These results are consistent with ours, when our upper plate is sufficiently weak (the trench follows the dynamics of the slab). However, ‘our’ trench motion remains constrained by that of the upper plate when it is sufficiently strong.

## 4.4 Comparison with the Andes and the Northwestern Pacific

The styles and cyclicity that were modelled here remain to be compared with specific subduction zones, measured plate velocities and tomographic images. Tomographic models indicate that on the Earth, slabs either stagnate in the mantle transition zone, or sink into the lower mantle (e.g. Fukao *et al.* 2001). However, because of poor resolution, results remain tricky to interpret in terms of slab geometry. Along the South American margin for instance, Li *et al.* (2008) showed that the slab crosses the 660 km depth transition zone in the Central Andes, and that it rather lies and even vanishes over this transition zone towards Southern Chile, coherently with its rejuvenation. In contrast, the other study by Fukao *et al.* (2009) (and Obayashi, personal communication, 2011) indicated that the Nazca slab remains above the 660 km discontinuity all along the Chilean

margin. Ongoing improvements to tomographic data are needed to determine slab geometries with respect to this transition zone.

Episodic magmatism along the Chilean Andes indicates that evolutionary cycles occurred during the Andean orogeny (Haschke *et al.* 2002). Slab steepening and shallowing with episodic durations of 30–40 Ma may explain eastward migration of the volcanic arc through time between 200 and 0 Ma, and periods of volcanic quiescence. Moreover, periods of shortening and extension in the South American continent are well identified, and include an Upper-Cretaceous ‘Peruvian’ shortening event that was widely observed from Ecuador to Patagonia around 80 Ma (e.g. Jaillard & Soler 1996; Ramos 2009), an Eocene–Oligocene (~39–49 Ma) ‘Incaic’ phase at different places along the Peru–Chile margin (Jaillard & Soler 1996) and the still active broad ‘Quechua’ phase from ~26 Ma (Pardo-Casas & Molnar 1987). Each of these phases are separated by volcanic gaps of ca. 5–10 Ma, which are thought to occur when the slab has a shallow dip and the mantle wedge is cold (Jaillard & Soler 1996; Haschke *et al.* 2002). These cyclic properties of Andean magmatism do not correlate with variations in plate motions, and had previously been interpreted as resulting from cyclic delamination of the overthickened continental crust (e.g. DeCelles *et al.* 2009; Ramos 2009).

Our modelled mode 1a cyclic behaviour offers a new explanation for the processes controlling the variations in dip of the Nazca slab. Along the Chilean margin, the present-day plate velocities are $v_{op} \simeq 4.3$ cm a⁻¹ and $v_{sp} \simeq 2.9$ cm a⁻¹ and a relative velocity equal to $v_s = v_{op} + v_{sp} \simeq 7.2$ cm a⁻¹, according to Nuvel-1A values (Pardo *et al.* 2002). Given these velocities and a standard plate viscosity $\eta = 1 \times 10^{24}$ Pa s, our numerical model produces cycles with a period of about 22 Ma and a shallow dip duration of ca. 4.3 Ma (Fig. 8). Episodes of a shallow slab dip at a minimum dip angle of 10° can explain the gaps in the observed volcanic activity. In our model the maximum dip angle reaches 48° (Fig. 8), whereas the steepest dip of the Nazca plate approaches 35° in Northern Chile. Our steep slab dips can be explained by the absence of a viscous asthenospheric mantle, and the lack of viscous drag forces acting on the slab. On the other hand, the shallow observed dip could indicate that the dip angle of the Nazca slab is not at its maximal possible value at present.

The period of modelled subduction cycles and the duration of a shallow slab are both a little shorter than those estimated for the Southern Andes (Haschke *et al.* 2002). The model assumptions of a constant slab viscosity, constant plate thickness and constant plate velocities through time are for an idealized system. A more accurate validation is beyond the reach of this study, since it requires to account also for climatic processes, 3-D effects, and other thermo–chemical–mechanical processes that have already been recognized to have a role in the Andean orogenesis. Despite these differences, our modelled values are within the range of the episodic variations estimated for the Andes. Haschke *et al.* (2002) argued that when the Nazca plate slab was steep in Central and Northern Chile, the South American margin underwent extension, and that a shallowing slab dip was concurrent with compression in the margin. Our numerical models show the same trend of strain in the overriding plate relating to slab dip. Therefore, cycle mode 1a might be a possible explanation for these observations. Further tomographic images of this part of the Earth’s upper mantle will help to determine whether slab folding is the primary mechanism that controls the slab dip and episodic deformation in this region.

The geometry of the western Pacific subduction zones has been well imaged with tomography models. The western Pacific slab lies extensively along the 660 km discontinuity underneath China, at
latitudes of 39–43°N (Fukao *et al.* 2009). The slab lying on the discontinuity appears to have a thickness twice that of the upper part towards the surface, which could be a result of slab folding and piling up. Naturally, such an assessment on the slab thickness must be taken with great caution, due to the uncertainties and trade-offs between spatial resolution and velocity contrasts in tomographic methods. According to recent kinematic estimations (DeMets *et al.* 2010), the overriding Eurasian Plate has an eastward velocity of $v_{op} \approx 20–24$ mm yr$^{-1}$, and the Pacific Plate subducts at a westward velocity of $v_{sp} \approx 63–69$ mm yr$^{-1}$ in this region. These values fall within the conditions of our cyclic mode 1a. Further investigation is required to accurately link slab folding cyclicity with any reported periodicity in the deformation of the Eurasian overriding plate.

## 5 CONCLUSION

We have developed a simplified 2-D numerical model that compares remarkably well with analogue models of slab deposition above the 660 km depth discontinuity (Heuret *et al.* 2007; Guillaume *et al.* 2009), and that shows similar styles of subduction. The lack of an asthenospheric mantle in our simulations leads to higher dip angles than for analogue models. However, this difference does not affect the occurrence of the four styles of subduction that are obtained in this study.

The style of subduction is controlled by the overriding plate velocity. The motion of the overriding plate towards the trench ($v_{op} > 0$) leads to style 1, with a slab that bends forwards as it lies on the 660 km depth discontinuity, whereas an overriding plate that is fixed or moves away from the trench leads to style 2, with a slab that bends backwards as it lies progressively on the discontinuity.

The occurrence of cyclicity within styles 1 and 2 is controlled by the relationship between both the upper plate and the slab velocities ($v_{op}$ and $v_{sp}$). When the absolute value of the overriding-plate velocity is equal to the subduction velocity ($|v_{op}| = v_s$), the process is steady-state. Otherwise, when $|v_{op}| < v_s$, cycles of slab folding occur. When $|v_{op}| > v_s$ the slab is stretched and thinned.

The occurrence of folds increases (folds become closer to each other) with a decreasing $v_{op}/v_s$ ratio at a given viscosity, or conversely when the viscosity of the subducting plate decreases at given ($v_{sp}, v_{op}$).

Changes in the stress regime of the overriding plate occur without modification of the boundary conditions, and are a consequence of folding cycles. When the dip of the slab shallows, the overriding plate undergoes compression and shows a positive change in topography. When the dip of the slab steepens, the overriding plate undergoes little compression (style 2) or extension (style 1), and a negative change in topography occurs.

The strength of the overriding plate exerts a first-order control on the trench motion and folding cyclicity, when all other conditions are fixed. The presence of a low-viscosity block in the upper plate acts as a buffer that impedes high stress storage and enhances cyclic compression and extension in the overriding plate, while the trench motion becomes more sensitive to slab behaviour alone. Therefore, the cyclicity of slab folding and trench motion both strongly depend on the strength of the overriding plate.

Numerical experiments with present-day velocities that correspond to the Central Andes margin evolve according to mode 1a (occurrence of forward folding cycles). A model with appropriate plate velocities produces a periodic deformation regime in the over-

riding plate that falls within the range of field observations (e.g. Kay & Mpodozis 2002; Haschke *et al.* 2002).

Modelling long-term subduction dynamics requires that the complexities of temperature and compositionally dependent rheologies are accounted for, including the solid-plate like behaviour at the top surface and decreasing viscosity towards the 660 km depth transition zone. It is also probably important to take into account surface processes as well as 3-D slab geometries, as they have been recognized to have an important role in the long-term building of the Andes (e.g. Lamb & Davis 2003).

## 6 ACKNOWLEDGMENTS

We thank B. Guillaume and J. Martinod for their helpful remarks and information about their analogue models. Many thanks to M. Obayashi and Y. Fukao for discussions and to G. Nolet for comments on this manuscript. J. Tetreault and an anonymous reviewer are thanked for their very constructive comments and helpful suggestions.

## REFERENCES

Arcay, D., Tric, E., Doin, M.-P., Bousquet, R. & de Capitani, C., 2006. Overriding plate thinning in subduction zones: localized convection induced by slab dehydration, *Geochem. Geophys. Geosyst.*, **7**, 1–26.

Bellahsen, N., Faccena, C. & Funiciello, F., 2005. Dynamics of subduction and plate motion in laboratory experiments: insights into the ‘plate tectonics’ behavior of the earth, *J. geophys. Res.*, **110**, doi:10.1029/2004JB002999.

Berger, A., Jouanne, F., Hassani, R. & Mugnier, J.-L., 2004. Modelling the spatial distribution of present-day deformation in Nepal: how cylindrical is the main Himalayan thrust in Nepal? *Geophys. J. Int.*, **156**, 94–114.

Billen, M., 2008. Modeling the dynamics of subducting slabs, *Annu. Rev. Earth Planet. Sci.*, **36**, 325–356.

Billen, M., 2010. Slab dynamics in the transition zone, *Phys. Earth planet. Inter.*, **183**, 296–308.

Bonnardot, M.-A., Hassani, R. & Tric, E., 2008a. Numerical modelling of lithosphere-asthenosphere interaction in a subduction zone, *Earth planet. Sci. Lett.*, **272**, 698–708.

Bonnardot, M.-A., Hassani, R., Tric, E., Ruellan, E. & Regnier, M., 2008b. Effect of margin curvature on plate deformation in a 3-D numerical model of subduction zones, *Geophys. J. Int.*, **173**, 1084–1094.

Buiter, S.J.H., Govers, R. & Wortel, M. J.R., 2001. A modelling study of vertical surface displacements at convergent plate margins, *Geophys. J. Int.*, **147**, 415–427.

Capitanio, F., Morra, G. & Goes, S., 2007. Dynamic models of downgoing plate buoyancy driven subduction: subduction motions and energy dissipation, *Earth planet. Sci. Lett.*, **262**, 284–297.

Capitanio, F., Stegman, D., Moresi, L. & Sharples, W., 2010. Upper plate controls on deep subduction, trench migrations and deformations at convergent margins, *Tectonophysics*, **483**, 80–92.

Chemenda, A., Lallemand, S. & Bokun, A., 2000. Strain partitioning and interplate friction in oblique subduction zone: constraints provided by experimental modeling, *J. geophys. Res.*, **105**, 5567–5581.

Chéry, J., Zoback, M. & Hassani, R., 2001. An integrated mechanical model of San Andreas fault in central and northern California, *J. geophys. Res.*, **106**(B10), 22 051–22 066.

Christensen, U.R., 1996. The influence of trench migration on slab penetration into the lower mantle, *Earth planet. Sci. Lett.*, **140**, 27–39.

Ciskova, H., Van Hunen, J., Van den Berg, A.P. & Vlaar, N.J., 2002. The influence of rheological weakening and yield stress on the interaction of slabs with the 670 km discontinuity, *Earth planet. Sci. Lett.*, **199**, 447–457.

Cundall, P. & Board, M., 1988. A microcomputer program for modeling large strain plasticity problems, in *Numerical Methods in Geomechanics*, pp. 2101–2108, ed. Swoboda, C., Balkema, Rotterdam.
DeCelles, P., Ducea, M., Kapp, P. & Zandt, G., 2009. Cyclicity in cordilleran orogenic systems, *Nat. Geosci.*, **2**, 251–257.

DeMets, C., Gordon, R. & Argus, D., 2010. Geologically current plate motions, *Geophys. J. Int.*, **181**, 1–80.

Duretz, T., Gerya, T. & May, D., 2011. Numerical modelling of spontaneous slab breakoff and subsequent topographic response, *Tectonophysics*, **502**, 244–256.

Fukao, Y., Widiantoro, S. & Obayashi, M., 2001. Stagnant slabs in the upper and lower mantle transition region, *Rev. Geophys.*, **39**, 291–323.

Fukao, Y., Obayashi, M., Nakakuki, T. & the Deep Slab Project Group, 2009. Stagnant slabs: a review, *Annu. Rev. Earth planet. Sci.*, **37**, 19–46.

Funiciello, F., Faccenna, C., Giardini, D. & Regenauer-Lieb, K., 2003a. Dynamics of retreating slabs: 2. Insights from three-dimensional laboratory experiments, *J. geophys. Res.*, **108**(B4), 2207, doi:10.1029/2001JB000896.

Funiciello, F., Morra, G., Regenauer-Lieb, K. & Giardini, D., 2003b. Dynamics of retreating slabs: 1. Insights from two-dimensional numerical experiments, *J. geophys. Res.*, **108**, 2207, doi:10.1029/2001JB000896.

Funiciello, F., Facenna, C. & Giardini, D., 2004. Role of lateral mantle flow in the evolution of subduction systems: insights from laboratory experiments, *Geophys. J. Int.*, **157**, 1393–1406.

Gerbault, M., 2000. At what stress level is the Indian ocean lithosphere buckling? *Earth planet. Sci. Lett.*, **178**, 165–181.

Gerbault, M., Cembrano, J., Mpodozis, C., Farias, M. & Pardo, M., 2009. Continental margin deformation along the Andean subduction zone: thermo-mechanical models, *Phys. Earth planet. Inter.*, **177**, 180–205.

Goes, S., Capitanio, F. & Morra, G., 2008. Evidence of lower-mantle slab penetration phases in plate motions, *Nature*, **451**, doi:10.1038/nature06691.

Goes, S., Capitanio, F., Morra, G., Seton, M. & Giardini, D., 2011. Signatures of downgoing plate-buoyancy driven subduction in cenozoic plate motions, *Phys. Earth planet. Inter.*, **184**, 1–13.

Got, J.-L., Monteiller, V., Monteux, J., Hassani, R. & Okubo, P., 2008. Deformation and rupture of the oceanic crust may control growth of Hawaiian volcanoes, *Nature*, **451**, 453–456.

Guillaume, B., Martinod, J. & Espurt, N., 2009. Variations of slab dip and overriding plate tectonics during subduction: insights from analogue modelling, *Tectonophysics*, **463**, 167–174.

Haschke, M., Scheuber, E., Gunther, A. & Reutter, K.-J., 2002. Evolutionary cycles during the Andean orogeny: repeated slab breakoff and flat subduction? *Terra Nova*, **14**, 49–55.

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**(B8), 17 951–17 965.

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., 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.

Hirth, G. & Kohlstedt, D., 2003. Rheology of the upper mantle and the mantle wedge: a view from the experimentalists, **138**, 93–108.

Huc, M., Hassani, R. & Chéry, J., 1998. Large earthquake nucleation associated with stress exchange between middle and upper crust, *Geophys. Res. Lett.*, **25**, 551–554.

Jaillard, E. & Soler, P., 1996. Cretaceous to early paleogenetectonic evolution of the northern central andes (0-18s) and its relations to geodynamics, *Tectonophysics*, **259**, 41–53.

Jarrard, R.D., 1986. Relations among subduction parameters, *Rev. Geophys.*, **24**(2), 217–284.

Kay, S. & Mpodozis, C., 2002. Magmastism as a probe to the Neogene shallowing of the Nazca plate beneath the modern Chilean flat-slab, *South Am. Earth Sci.*, **15**, 39–57.

King, S.D., 2001. Subduction zones: observation and geodynamic models, *Phys. Earth planet. Inter.*, **127**, 9–24.

Lallemand, S., Heuret, A., Faccenna, C. & Funiciello, F., 2008. Subduction dynamics revealed by trench migration, *Tectonics*, **27**, doi:10.1029/2007TC002212.

Lamb, S. & Davis, P., 2003. Cenozoic climate change as a possible cause for the rise of the Andes, *Nature*, **425**, 792–797.

Lesne, O., Calais, E., Deverchere, J., Chéry, J. & Hassani, R., 2000. Dynamics of intracontinental extension in the north Baikal rift from two-dimensional numerical deformation modeling, *J. geophys. Res.*, **105**, 21 727–21 744.

Li, C., der Hilst, R.D.V., Engdahl, E. & Burdick, S., 2008. A new global model for p wavespeed variations in Earth's mantle, *Geochem. Geophys. Geosyst.*, **9**, doi:10.1029/2007GC001805.

Neves, S.P., Tommasi, A., Vauchez, A. & Hassani, R., 2008. Intraplate continental deformation: influence of a heat-producing layer in the lithospheric mantle, *Earth planet. Sci. Lett.*, **274**, 392–400.

Ord, A. & Hobbs, B., 1989. The strength of the continental crust, detachment zones and the development of plastic instabilities, *Tectonophysics*, **158**, 269–289.

Pardo, M., Comte, D. & Monfret, T., 2002. Seismotectonic and stress distribution in the central Chile subduction zone, *J. South Am. Earth Sci.*, **15**, 11–22.

Pardo Casas, F. & Molnar, P., 1987. Relative motion of the Nazca (Farallon) and South American plates since late cretaceous time, *Tectonics*, **6**, 233–287.

Poliakov, A.N. & Podladchikov, Y., 1992. Diapirism and topography, *Geo-phys. J. Int.*, **109**, 553–564.

Ramos, V., 2009. Anatomy and global context of the Andes: main geologic features and the Andean orogenic cycle, in *Backbone of the Americas: Shallow Subduction, Plateau Uplift, and Ridge and Terrane Collision*, Geological Society Memoirs Vol. 204, pp. 31–65, eds Kay, S., Ramos, V. & Dickinson, W., Geological Society of America, Boulder, CO.

Ribe, N., 2010. Bending mechanics and mode selection in free subduction: a thin-sheet analysis, *Geophys. J. Int.*, **180**, 559–576.

Royden, L. & Husson, L., 2006. Trench motion, slab geometry and viscous stresses in subduction systems, *Geophys. J. Int.*, **167**, 881–905.

Schellart, W., 2005. Influence of the subducting plate velocity on the geometry of the slab and migration of the subducting hinge, *Earth planet. Sci. Lett.*, **231**, 197–219.

Schellart, W., 2008. Subduction zone trench migration: slab driven or overriding-plate-driven? *Phys. Earth planet. Inter.*, **170**, 73–88.

Schellart, W., 2009. Evolution of the slab bending radius and the bending dissipation in three-dimensional subduction models with a variable slab to upper mantle viscosity ratio, *Earth planet. Sci. Lett.*, **288**, 309–319.

Schellart, W.P., 2010. Introduction to the special issue on convergent plate margin dynamics, *Tectonophysics*, **483**, 1–3.

Stegman, D.R., Farrington, R., Capitanio, F.A. & Schellart, W.P., 2010a. A regime diagram for subduction styles from 3-D numerical models of free subduction, *Tectonophysics*, **483**, 29–45.

Stegman, D.R., Schellart, W.P. & Freeman, J., 2010b. Competing influences of plate width and far-field boundary conditions on trench migration and morphology of subducted slabs in the upper mantle, *Tectonophysics*, **483**, 46–57.

Underwood, P., 1983. Dynamic relaxation, in *Computational Methods for Transient Analysis*, pp. 245–265, eds Belytschko, T. & Hughes, T., Elsevier Science, New York, NY.

Van Hunen, J., Funiciello, F. & Faccenna, C., 2008. Slab stiffness control of trench motion: insight from numerical models, *Geochem. Geophys. Geosyst.*, **9**, doi:10.1029/2007GC001776.

Vanbrabant, Y., Jongmans, D., Hassani, R. & Bellino, D., 1997. An application of two-dimensional finite-element modelling for studying the deformation of the Variscan fold-and-thrust belt (Belgium), *Tectonophysics*, **309**, 141–159.

Yamato, P., Husson, L., Braun, J., Loiselet, C. & Thieulot, C., 2009. Influence of surrounding plates on 3d subduction dynamics, *Geophys. Res. Lett.*, **36**, doi:10.1029/2008GL036942.

Yoshioka, S. & Naganoda, A., 2010. Effect of trench migration on fall of stagnant slabs into the mantle, *Phys. Earth planet. Inter.*, **183**, 321–329.